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CHAPTER I 


INTRODUCTION 


Instability of viscous flows along concave surfaces, like that 
between rotating concentric cylinders, is controlled by the balance 
of induced centrifugal forces and other forces acting on the flow. 

The criteria for instability for such systems was first postulated by 
Rayleigh (1916) for inviscid flows, and examined both theoretically and 
experimentally by Taylor (1923) for circular Couette flow, while Goertler 
(1954) investigated the mode of motion in a boundary layer along a con- 
cave surface. Centrifugal forces induced by curvature effects in these 
systems lead to the instability of the flow in the form of counter- 
rotating, vortex-like disturbances. 

Earlier transition measurements in the incompressible boundary 
layer next to a concave surface, indicate a steady three-dimensional 
vortex-like disturbance with a spanwise periodicity which develops 
according to the linearized theory. Gregory and Walker (1956) were the 
first to observe traces of these vortices by using the China-Clay tech- 
nique, followed by Aihara (1962) and Tani and Sakagami (1962) using 
colored liquids and smoke threads, and Aihara (1961) and Tani (1961) 
using hot-wire measurements. Wortmann (1964) used the telurium method 
to visualize these vortices in a water tunnel. Bippes (1978) and 
Bippes and Goertler (1972) presented detailed observations of these 
vortices using the hydrogen-bubble technique. At compressible speeds, 
evidence of the vortex-like disturbances has been observed by Persen 
(1968) and Ginoux (1970) in quasi two-dimensional flows in regions of 
separated flow reattachments. Zakkay and Calarese (1972) observed 
the presence of these vortices in a hypersonic turbulent 
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boundary layer over an axi symmetric configuration with adverse pressure 
gradient. In their experiments in two Mach 5 nozzles, Beckwith and 
Holley (1981) showed by using oil flow patterns that these vortices 
persisted to the nozzle exit and that the vortices were involved in the 
transition process. 

Available experimental evidence has shown that the counter-rotating 
vortices affect indirectly the transition from laminar to turbulent flow 
on a concave surface. As steady vortices, they may not lead to transi- 
tion by themselves, but when their growth becomes strong enough, they do 
cause transition by an unknown mechanism. Transition over a concave sur- 
face occurs at Reynolds numbers that are lower than those for flow over 
a convex or flat surfaces (Clauser and Clauser, 1937). Moreover, the 
Goertler number (the nondimensional parameter of interest in this type 
of problems) has to reach a certain critical value before transition 
takes place over a constant curvature surface (Liepmann, 1945). Tarn* 
indicated that, although their growth is small these vortices cause 
spanwise variation in the velocity field, thus modifying the development 
of unstable waves. Wortmann (1969) concluded that transition was lead 
by a secondary instability that he observed following the appearance of 
the counter-rotating vortices. On the other hand, the observations of 
Bippes (1978) and Aihara (1976) of a sinusoidal motion of the vortex 
axes before turbulence sets in lead Aihara to correlate that with a non- 
linear theory. Nayfeh (1979) studied the effect of these vortices on 
Tollmien-Schlichting waves, and showed that these vortices have a strong 
tendency to amplify three-dimensional waves having a spanwise wavelength 
that is twice the wavelength of the vortices. 

The instability of boundary layer flow on curved surfaces was first 



demonstrated theoretically by Goertler (1954) for incompressible flows, 
where he showed that a system of counter-rotating vortices were formed 
parallel and oriented in the streamwise direction. These vortices are 
referred to as Taylor-Goertler vortices or also as Goertler vortices. 

In his analysis, Goertler assumed the boundary layer to be parallel, 
streamwise curvature to be constant with the distance normal to the 
surface, and the vortices to be confined to the boundary layer. In an 
attempt to relax these assumptions, several investigators followed 
Goertler and extended his analysis for incompressible flows. A detailed 
review of these efforts is given by Herbert (1976) and Floryan and Saric 
(1979). 

Various theoretical investigations were performed to provide an 
accurate mathematical model of the instability nearer to the physical 
reality. Various investigators used the body oriented coordinate sys- 
tem 1 with some assumptions regarding the variation of curvature in 
the direction normal to the flow. The governing equations written, in 
the body oriented coordinate system contain the terms (1 - Ky)“ n , where K 
is the curvature of the surface, and y is the coordinate normal to the 
curved surface. These terms present a singularity when Ky = 0(1). Goertl 
(1954) assumed that these terms can be replaced by 1, that is stream- 
line curvature is constant at any normal distance. Smith (1955) and 
Kahawita and Meroney (1977) expanded these terms binomial ly and kept 
the first two terms, i.e. (1 + nKy), thus effectively transferring the 
singularity to infinity. Floryan and Saric (1979) gave a brief descrip- 
tion of different approaches used in the past using body oriented coor- 
dinate system and utilized a new system of coordinate based on the stream- 
lines and potential lines of the inviscid flow over the curved surface. 



A typical illustration of this coordinate system is shown in Fig. 2 for 
a circular arc. The shape of the streamline is directly related to the 
system of coordinates. Consequently, different treatments of the stream- 
line curvature are expected to influence the stability characteristics 
of the flow because they are equivalent to changes in the outer flow 
conditions. Larger rates of decay of curvature outside the boundary 
layer causes a reduction of the driving centrifugal forces in the outer 
flow. A similar situation exists by decreasing the streamwise extent of 
the curved flow region (see Herbert, 1976). Both factors considerably 
stabilize the flow. 

The effect of boundary-1 ayer growth was introduced by several 
researchers (e.g. Smith, 1955; Kahawita and Meroney, 1977; Herbert, 

1976; Floryan and Saric, 1979; Ragab and Nayfeh, 1980). The inclusion 
of the normal velocity component of the basic flow in the analysis dras- 
tically changes the location of the neutral curve. Moreover, the stream- 
wise variation of the normal velocity component which appears in the 
leading order stability equations has a strong influence on the stabil- 
ity characteristics. While Smith incorporated the normal velocity com- 
ponent of the basic flow with some of the higher order curvature 
terms, Floryan and Saric based the scaling of the normal disturbance 
velocity component on a viscous scale that contributed a leading order 
effect of the normal velocity component of the basic flow and its stream- 
wise variation on the stability analysis. 

In spite of the extensive investigations of Goertler instability in 
incompressible flows, there have been only a limited number of studies 
on the effect of compressibility of the basic flow on this type of distur- 
bance. The compressible linear stability theory now available (Hammer! in, 
Aihara, 1961; and Kobayashi and Kohama, 1977) for the development of these 



vortices Is not reliable and far from being complete. It 5 

neglects the normal velocity component of the basic flow that proved 
to have a profound effect in reducing the stability of incompressible 
flows. Existing compressible theories treat only the neutral stability 
case which is of limited importance regarding vortex development, pos- 
sible nonlinear interactions, and transition correlation. 

In his compressible stability analysis of boundary layer flow along 
a concave surface, Hammerlin expanded the disturbance equations in power 
series of the square of the freestream Mach number. He kept terms to 
0(M«), and obtained a solution to the disturbance equations which is 
valid only for Moo « 1. He used a power law for the viscosity-temperature 
relation and kept Prandtl number constant in the calculations. His 
results Indicated that minimum critical Goertler number increases as 
Mach number increases. Aihara reduced the perturbation equations into 
two extreme cases of the freestream Mach number, namely, M^ « 1, and 
» 1. In his analysis, he assumed constant viscosity and Prandtl 
number. His results shows that minimum critical Goertler number de- 
creases as Mach number increases, that is, compressibility effect on 
the stability of the boundary layer is opposite to that predicted by 
Hammerlin. Kobayashi and Kohama treated the problem over a wide range 
of Mach numbers. They used Sutherland's formula for the temperature 
dependence of the viscosity and assumed constant Prandtl number and 
specific heat. They treated only the neutrally stable disturbances 
showing that the boundary layer becomes more stable (Goertler number 
increases) as Mach number increases. 

Kobayashi (1972,73,74) was the first to examine the effects of suc- 
tion on the stability characteristics of a laminar incompressible boun- 
dary layer along curved surfaces. Although his analysis excluded the 
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effect of normal component of mean velocity due to boundary layer 
growth, it gave insight into the effect of the presence of this velocity 
component due to suction. Using homogeneous suctions Kobayashi found 
that the laminar boundary layer is stabilized (critical Goertler number 
increases). However, the changes in the critical Goertler number re- 
mains much less than the changes in the critical Reynolds number for 
the Tollmien-Schlichting instabilities due to homogeneous suction 
(Hughes and Reid, 1965). So, Goertler instability will probably pre- 
dominate in the laminar boundary layer along the concave wall with suc- 
tion. Floryan and Saric (1979) and Floryan (1980) included the normal 
component of the mean velocity in their analysis and examined the case 
of self similar suction and came to the same conclusion. 

DiPrima and Dunn (1956) examined the effect of cooling and heating 
on the stability characteristics of laminar liquid boundary layer over 
curved surface. Since their analysis was done for liquids, they neglec- 
ted the viscous dissipation and variations in density and came to a 
conclusion that heating or cooling has very slight influence on Goertler 
instability. For compressible boundary layer, Kobayashi and Kohama 
(1977) found that for isothermal walls the ratios of wall to freestream 
temperature has less effect on the critical value of Goertler number as 
Mach number is increased. 

In this article, a basic approximation to a compressible linear 

stability theory is developed for three-dimensional longitudinal type vortices 

in two-dimensional compressible boundary layers along curved surfaces. 

The effect of compressibility on the critical stability limit, growth 
rates, and amplitude ratios of the vortices is evaluated over a range 
of Mach numbers from 0 to 5. The effect of boundary layer growth is 



included in the analysis. The effect of wall cooling and suction on 
the development of Goertler velocities in a compressible boundary layer 
is examined. In Chapter II, a formulation of the stability problem is 
introduced. In Chapter III, the method of solution and numerical pro- 
cedures are outlined. A discussion of the numerical results is given 
in Chapter IV and conclusions are in Chapter V. 
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CHAPTER II 


PROBLEM FORMULATION 


The Spatial three-dimensional stability of laminar incom- 
pressible and compressible two-dimensional boundary layers along 

a slightly curved wall is considered. The wall curvature 
is in the direction of the flow and its variation is assumed to be 
weak to avoid a nonuniform Mach number distribution along the 
wall due to the presence of shock waves that might occur in the case 
of rapid changes. 

The flow field is governed by the Navier-Stokes, energy, 
continuity and state equations written in an orthogonal 
curvilinear coordinate system. The local curvature of the 
streamlines enters the field equations through the appropriate 
coordinate system. Following Floryan and Saric and Ragab and Nayfeh 

a coordinate system (x, y, z) based on the streamlines and potential 
lines of the inviscid flow over a curved surface is used in this 
analysis. Here x and y are in the direction of stream lines and po- 
tential lines, respectively, and z is the corrdinate normal to the 
y plane. Figure 3 shows the direction of the development of 
Goertler vortices in this coordinate system. This coordinate 
system has the advantage of a body oriented coordinate 
system in the wall region and decays to a rectangular Cartesian 
coordinate system away from the wall. 

The compressible field equations written in this coordinate 
system are presented in Appendix A. The local curvature of the 
streamlines enters the field equations through the metric 
coefficients. They are determined from the definition of the 



arc length to be h , h and 1 in the x, y and z directions, respectively . 
The field equations given in Appendix A are used to formulate the 
disturbance equations. 

We consider the basic state to be two-dimensional, viscous, 
compressible flow over a slightly curved surface. The field 
equations are made dimensionless using a reference length 

L* = (v* x*/U*) (1) 

and U* as reference velocity and(p*u* 2 )oo as reference 
pressure, where * indicates a dimensional quantity, and v is 
the kinematic viscosity of the fluid. The thermodynamic and 
transport properties of the fluid are made dimensionless using 
their corresponding free stream values. With these definitions 
the characteristic Reynolds number becomes 

R = (USL*/>&) ( 2 ) 

Vie define a small viscous parameter e , and a small curvature 
parameter k as 

e = 1/R, k =(L*K*)l/2 


where K* is the curvature of the wall. The metric coefficient 
h is related to the curvature parameter k by 


k 2 =-h /h 2 at y = 0 


( 4 ) 


2.1 The Mean Flow 


Van Dyke (1960,62) using a body oriented coordinate system, 
showed for incompressible and compressible flows that the leading 
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order approximation of the boundary layer equations is the familiar 
boundary layer equations on a flat plate. 

In terms of the small parameters e and k, Floryan (1980) showed 
that for incompressible flows the leading order approximation of the 
boundary layer equations (terms of order e° and k°), when e and k are 
of the same order, is the conventional boundary layer equations over 
a flat plate. Modification of the boundary layer profiles due to the 
curvature comes mainly through the normal pressure gradient. This 
pressure gradient is 0(k) and it will enter higher order boundary layer. 

Because the aim of this study is to provide a basic approximation 
for the stability of a compressible boundary layer over a curved sur- 
face, therefore the basic approximation for the mean flow is required. 
The compressible boundary layer flow over a flat plate is considered 
to provide a basic approximation for the stability of the compressible 
boundary layer over a curved surface. The mean flow profiles are cal- 
culated for an adiabatic or isothermal flat plate. In case of suction, 
a similar suction parameter y is defined as 

(pV) w (5) 

y - R 

(p*u*)« 

and is introduced to the boundary layer equations (see Appendix B). 

Here w indicates wall conditions, and R is given in Eq. (2). The 
fluid is considered to be a perfect gas with all the thermodynamic 
and transport properties function of temperature. The mean 
viscosity is related to mean temperature through Sutherland for- 
mula. The flow stagnation temperature is kept constant and equal 
to 310K for all Mach numbers under study. 



2.2 The Disturbance Flow 


A steady three-dimensional small disturbance is superposed on 
each mean -flow quantity in the following form 


u(x,y ,z) = U(x,y) + u(x,y,z) (6a) 

v(x,y ,z) = V(x,y) + v(x,y,z) (6b) 

w(x,y,z) - 0 + w(x,y,z) (6c) 

p(x*y,z) = P(x,y) + p(x,y,z) (6d) 

6(x,y,z) = 0(x,y) +e(x,y,z) (6e) 

p(x,y,z) = p(x,y) + p(x,y,z) (6f) 

p(x,y,z) = p(0) + ^ 0(x,y,z) (g g ) 


In the above equations the order of magnitude of the normal mean veloc- 
ity V is smaller than U by the Reynolds number R, and the disturbance 
quantities are made dimensionless using 

u = u*/l£» v = v*/RU*, w = w*/RU*. 

P = P*/R 2 p*U* 2 , e = 0V0*. p * p*/p* ^ 

The order of magnitude of the disturbance quantities v and w differ 
from the disturbance U by the Reynolds number R. The need for differ- 
ent disturbance velocity scaling was recognized by DiPrima and Stuart 
and observed by Bippes, Bippes and Goertler.and Wortmann in 
their experiments. Floryan and Saric and Ragab and Nayfeh used this 
type of scaling and as a result the normal velocity component of the 
mean flow affected the leading order stability analysis. With this 
scaling the disturbance motion varies in the 
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x coordinate according to a new scale. 

x = x*/RL* (8) 

Boundary- layer theory combined with the assumed scales, gives the x vari 
atl&n' of mean-flow quantities in terms of the new scale. 

Substituting Eq. (6) into the dimensionless field equations, sub- 
tracting the mean flow, linearizing the equations, and keeping the lead- 
ing order terms, that is terms 0(1), the following disturbance equa- 
tions are obtained: 


x-momentum 


y-momentum 


z + Iv (pU y ) y + 5 u y v 
- [ ^- (uu x + vu y ) + (qi ) y )^ 0 - C u y e y = ° 

l(V x + 2UG Z )u -c My u x - (c + l)uu yx - V y 
+ |< Vv y^ yV zz- (c+ 2 )(vV yV + Py 

-[•^ UV x + VV y + u2g2) + (c + 1)iiU yx + cii y u x 
+ (C + 2){i3V y ) y + u x ul e - W y e x -[cyU x 
+ (c + 2)CvJ e y - cy y w z - (c + 1 )uw y2 = 0 


(9) 


( 10 ) 


z-momentum 


w x u z + + 1 ^“xz + M y v z + ( c + 1 ^ v yz " p z 

+ cy(U x + V y )e 2 - 1 w x + (c + 2)ww 2Z - 1 w y 


(ID 


+ (ww y ) y = 0 
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Energy 


I 0 X U - 2(Y- 1 )(&U yUy 4 V ‘[^ (US x + V V 
+ ( Y - lJM^tu ) 2 + ]KD0 y Ge + | e x - r 0 Z2 

J (13) 

+ (g-f P 0 y) e y ' f^Vy ’ 0 


By keeping the leading order terms in the disturbance equations, 
the parameters e and k appear only implicitly in the so called Goertler 
number G and consequently Eqs. (9) - (13) are configuration independent 
at this level of approximation. The parameters e and k will appear 
explicitly in a higher order stability analysis. The Goertler number 
G is defined by 

G » Rk (14) 

The other parameters that appear in. Eqs. (9) - (13) are the free stream 
Mach number M», Prandtl number r, and p = dp/de. The density disturbance 
is eliminated by using the equation of state. 
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Equations (9) - (13) constitute a system of homogeneous linear partial 
differential equations with variable coefficients. The type of distur- 
bance considered in this analysis is periodic in the z-coordinate, 
thus the dependence on z can be eliminated from the equations. The 
scale x may also be separated from the equations by considering the 
mean flow to be quasi parallel. Then the form of the disturbance and 
flow characteristics suggest the normal mode approach. The disturbance 


variables can be written as 

u = u(y) cos (8z) exp (/adx) (15a) 
v = v(y) cos (0z) exp (/adx) (15b) 
w - w(y) sin (0z) exp (/adx) (15c) 
p = p(y) cos (0z) exp (/adx) (15d) 
0 = 6(y) cos (0z) exp (/adx) (15e) 


where 8 is the dimensionless wave number in z-coordinate, and a is the 
spatial growth rate. The variables u, v, w, p, and 6 should not be con- 
fused with the total flow variables defined in Eq. (6) which will not be 
referred to later on. 

Substituting Eq. (15) into Eqs. (9) - (13) we obtain 

£i{aU + U x )+ve 2 Ju+ (|-p y )u y - vO yy + |u y v 

-[4( u Ux + VU y ) + (iiU y )^0-uU y 0 y = 0 ( 16 ) 


x-momentum 
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y-momentum 


z-momentum 


Continuity 


Energy 


^(V x + 2UG 2 ) - cay^J u - £(c + 1 )ay + y 1 u y 
+ |(v y + aU + ne Z )v (c + 2)yl V y 
- (c + 2)yv yy - cu y Bw - (c + 1 )y8w y + p y 


- -7(uv x + vv y + u 2 g 2 ) + (c + i )qu xy + cq y u x 

— 0 — * _ 


+ (c + 2)(pV y ) y + M x U y + auU y 


+ ( c + 2 ) jjV y 0 y = 0 


«4 — 

0- cyll 

U 


(17) 


3 j^( c + 1 )op + u + Bp y v + (c +. 1 )Buv y + |^(c + 2 ) 3 2 m 

+ - (u y - |)w y - yw yy - $p + ceq(U x 

+ V y )e = o 


(18) 


l/ a _ 1 0 )u_J_Q y + lj + £ w + — (— 0 + — 0 U9) 

Q^° 0 u x' u q 2 u y v + 0 v y + O w + 0 2 0 u x 0 u y 


A U /S 

A * A 


- W^ 0 -^ 6 * = 0 


l0 x G - 2 (Y-l)HfpU y G y + l0 y v-[^-(U0 x + V0 y ) (20) 

+ (Y-i)Mfp(u y ) 2 + i(qo y ) y - 
+ r^®y " r'V^y * T^yy = 0 
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In the limiting case of Moo = 0 the above equations reduce to the 
incompressible equations developed by Floryan and Saric, and reduce to the 
basic system of equations developed by Ragab and Nayfeh. The above 
equations also reduce to the compressible equations of Kobayashi and 
Kohama for the neutral stability case (a = 0) and when all normal velo- 
city and streamwise variations of the mean flow are disregarded. 

Eq. (16) - (20) are supplemented with the appropriate boundary 
conditions, which are 

u = v - w = 0, and 

0 y + B6 = 0 at y = 0 (21) 

u,v,w,§ -*■ Oas y (22) 

The function B in the condition (21) depends on the thermal pro- 
perties of both the solid surface and the fluid, as well as the wave 
length of the disturbance. For an adiabatic surface the case where B -*■ 

A 

0(6y = 0 at y = 0 that is to consider "adiabatic disturbance") is con- 
sidered. However, for cases involving surface cooling the case where 

.A. 

B 00 (0 = 0 at y = 0 ) is considered. The effect of using the thermal 
boundary condition 0y(O) = 0 on the stability characteristics of the 
boundary layer over an adiabatic surface is discussed in section 4.3. 
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CHAPTER III 

COMPUTATIONAL AND NUMERICAL PROCEDURE 


Equations (16 - 22) constitute an eighth -order system of homoge- 
nous linear differential equations with homogeneous boundary conditions 
This forms an eigenvalue problem for the parameters 3, a, and G. 
Equations (16 - 22) are reduced to eight first-order equations in the 
form 


(Zm)y - 

8 

n -l a mn^n = 0 

m = 1 , 2, . . . 

8 

(23) 

h - Z3 

- Z 5 = Z 7 = 0 

at y = O,for 

isothermal wall 

(24) 

or Zj = Z 3 

It 

M 

03 

ii 

fNl 

ii 

o 

at y - O,for 

adiabatic wall 

(25) 

z l» z 3» 

z 5 , z 7 -0 

as y + » 


(26) 


where Z's are defined as 


Z 


1 


Z 


5 


u> “ u y* ^3 = v » ^4 “ P 
0, Zg = 0y, 2-i = w, Zq = w y 


(27) 


and a mp is a variable coefficient matrix, whose nonzero elements are 


a 12 “ 

1 



*21 * 

*>< U x 

+ aU) 

+ e 2 

a 22 = 

1(V _ 

y v 0 

M y ) 


a 23 = 

. 1 u 
u0 y 



a 25 * 

--f-L(uu + 
w Le 2 ' x 

v V ♦ eu yy * q y u y ] 
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a,, = -- u„ 

26 u y 

a 31 = i 0 x " a 


a 33 = 0 0 y 


l 35 “ t + K + V y> - |-( U0 x + V0 y> 


36 " 0 


a 37 “ 3 


41 




7ET I 2UG + 2ap + V x + | © x - aV 


rv 

- (c * 2)(- ay + 0 x y y + / 0 X + 0 xy n) 


a 42 = - oy + 


43 


-it 


.c + 2 r 

v— [ ( 


0 X - 2(y-l)M>U y 


aU + liOg + V y + | 0 y - (c + 2)(0 y y y 


+ U0 + 

M yy 


7T- © ) 

G y' 


l 45 


U 2 G 2 , c + 2 


0 


2“ + ^ aU y - | e y> + ^ 


+ (c + 2)[(MV y ) y -|(M0y) y ] + J(UVx + W y ) 

i[f< V0 y 

' o| j)[|- < c + 2 )( Vg 2 M0 y 

- ^ - Q l n [|( Ve y + U 0 X> - OU+ (Y- l)M 2 00 (U y ) 


+ (c + 1 + C x U y + cd y U x + 


+ U0 x> - U x - V y 


♦f e 2 


2V0 + 2U0 + 30 V 

yy xy Ju y v y 


m 

+ 2Q x U y - f( Q y) 2 • f 0 y 0 x + 0yU x - eu X y - 0V yy J 
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l 46 


rv 
y ' 0 u 


■l^[oU-fe y -^ VV 2v 
*] - £ ♦ * <*«, 


-Ip -2L 

v y v 

+ (c + 2)ijV t 


it 


a 47 = - 23 m + h |3V - (c + 2)Bu0, 


a 48 = PlJ 


a 56 = 1 


a 61 p0 ®x 


a 62 = 2 (y - l)M‘ru 


] 


a C o = -pr 0 

63 y0 y 


*»s ' <* - >,>, * J [f - jf<»„ * »V 

] 


- (y - i)M 2 p(u y ) 2 


l rv 

a 66 = p(~0 ‘ 'V " 
a 78 = 1 


a 81 = g( K + £ ^ J ' 0 x ) 


a 83 " 6 C y y + C 0 1 0 y> 
a 84 = 'p 


'85 


■**4 


oU + U x + V y * |< V9 y + U9 x> 


] 


+ £Mi(U + y ) 
p ' x y' 


a 86 ■ Hr- 6V 


a = ft 2 . oU 
a 87 6 p0 


. . I/V , 

a 88 ' p l 0 ■ ^y ' 
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The reduction of the governing equations to a system of eight first- 
order equations requires only the second derivatives of the basic-flow 
quantities. In the work of Kobayashi and Kohama the disturbance 
equations for neutral stability were reduced to three coupled equations, 
the first of second order in u,the second of fourth order in v,and the 


third of second order in 0»and then solved numerically by finite differ- 
ence. The coefficients in their equations contain fourth-order deriva- 
tives of the basic flow, which requires a higher accuracy of the basic 
flow solution than the method used here. 

3.1 Eigenvalues and Eigenvectors 

Outside the boundary layer, i.e. at y ^y e (e indicates the edge of 

the boundary layer), the coefficients in equations (23) become constant 

since u 3 v s 8 8 1, P s 3 constant, V = V e = constant, r = r e = 

constant and the x and y derivatives of all mean-flow quantities are 

zero except V x . The nonzero elements of the constant coefficient matrix 
* 




are: 

<2 = 1 

a*i “ a + 6 : 

a* 2 = V e 

* 

a 3 1 = “o 

a 3 s - o 

•« - V e 

* 

a 37 = “6 

★ 


«« ■ [2G 2 + V xe - oV e ] 


* 

a 42 ** -a 


at 3 .* -(a + e*) 
at 5 = G 2 + V xe - oV e 


* (C + 2)reVe ^B 2 * Oj 


a* 6 

= (c + 

2 )(c 

F + 

a 47 

= 8 V e 



* 

a 48 

=-$ 



a te 

= 1 



* 

a 6 5 

* B 2 + 

r e a 



r e v e 




* 1 



a* 
a 8 4 

* -0 



■ * 

a 8 5 

■ (c + 

1 ) 

00 

* 

a 6 6 

= (c + 

1 ) 

0V ( 

* ■ 
a® 7 

= 3 2 + 

a 


a 88 

v e 




e*e 
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Equation (23) with the constant coefficient matrix a mn has a solution 

that can be expressed in the form 
8 

Zm = E D mn c n exp(X n y) m = 1 , 2 ,. . . , 8 , y=y e ( 28 ) 

n=l 

where D mn is the characteristic eigenvector matrix, c^is arbitrary con- 
stant vector ,and X n are the eigenvalues of the matrix a mn . They are: 


A 2 = A 3 = 4 { - V e + t V e + 4 (° + 6 2 )] 1/2 } 

A 4 = 4 { ' r e V e + [r e V e + 4 (° r e + 6 2 )1 V2 > 

6 6 ® ® (29) 

A 5 - 6 

A6 = X 7 = -|{V e + (V 2 + 4(o + 6 2 )] 1/2 ) 
a 8 * 5 -(r e v e + [r 2 v 2 + 4(ar e + e 2 )l 1/2 J 

We consider the case in which \ lt X 2 , X 3 and X 4 are negative real numbers 
whereas the other four are positive real numbers. Only negative signs 
satisfy the boundary conditions (26). 

It should be noted that the eigenvalues \ 2 anc j a 6 are repeated. 

Also in the special case when a = 3V e four eigenvalues are repeated ,i .e. 


Ai — .X 2 = X 3 - X 4 - -3 (30) 

A 5 = Xg = X 7 = Xe = 3 

Evaluation of the eigenvectors requires special care when the multipli- 
city of the roots is equal to or greater than two. Eigenvectors cor' 
responding to the eigenvalues (29) are given In Appendix C and those cor- 
responding to the eigenvalues (30) are given in Appendix D. Eigenvectors are 
checked by reducing the constant coefficient matrix a* to a Jordan canon- 
ical form using the similarity transformation 

J = D _1 a*D (31) 

(where D" 1 is the inverse of the eigenvector matrix D). 

When A n are distinct (which is not the case here), the resulting matrix 
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J is a diagonal matrix with the diagonal elements as eigenvalues. In 
case of repeated roots, elements of the super diagonals corresponding 
to the repeated roots are nonzeros. 

3 . 2 Boundary Conditions 

The asymptotic boundary condition (26) as y-*® demand that the con- 
stants c 5 - c 8 be zero in the solution (28). This can be expressed as 
8 

£ 0’/, Z m = 0 m = 5, 6 , 7, 8 y = y e (32) 

n=l 

The boundary condition at the wall (24) or (25) can be written in the form 
8 

£ e mn Z n = 0 m = I-,- 2, 3, 4, y = 0 (33) 

n=l 

where the elements of the 4x8 matrix e mn are all zeros except eu = 

623 “ 6 36 = e 47 = 1 in case of adiabatic wall, and eu = e 2 3 = e 35 = 
e 47 = 1 in case of isothermal wall. 

3.3 Numerical Procedures 

Equations (23) with boundary conditions (32) and (33) form a two 
point boundary value problem which is solved numerically. We assign values 
for two of the parameters B, a and G, guess the third one, and use the 
boundary condition (32) to construct a linear combination of the general 
solution (28). A variable stepsize integrator, written by Scott and Watts 
(1977), based on the Runge Kutta-Fehlburg fifth-order formulas, coupled 
with orthonormalization is used to integrate equations (23) from y = y e 
to the wall. At the wall, the values of the independent solution vectors 
are linearly combined to satisfy all but one of the wall boundary conditions. 
The last wall boundary condition can only be satisfied by this combined 
solution when the exact 
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eigenvalue has been found. A Newton-Raphson procedure is used to de- 
termine this eigenvalue. In the results presented here, a value is, 
always assigned to o to locate contours of constant growth rate 
use either 3 or G as an eigenvalue and determine it to 0(10"^). A 
FORTRAN computer program for the stability analysis is given in Appendix F. 

3.4 Adjoint Problem 

The system of equations adjoint to (23) are formed to check the 
eigenvalue. The adjoint system of equations can be written as 
8 

(^m)y + a nm Z n = 0 m = 1, 2, ...,8 (34) 

n=l 

ZJ = Z4 = Z$ = Zg = 0 at y = 0, for adiabatic wall (35) 

or z £ = Z| = Zg = Zg = 0 at y = 0, for isothermal wall (36) 

Z 2 ’ Z 4 » Z 6 * Z 8 0 as y (37) 

Equations (34) - (37) are solved using the same procedures used to 

solve the regular equations (23) - (26). The eigenvectors of the 
adjoint problem are given in Appendix E. 

Eigenvalues calculated by Equations (23) and its adjoint (34) were 
checked, and found to agree to 0(1CT 4 ), and both were found to be in- 
dependent of the value of y e . However, the value of y e used for the 
adjoint problem (34) to achieve the required accuracy was always higher 
than the value of y fi used for the regular problem (23). It was also 
necessary to increase y e as Mach number increases. For example for the 
regular problem y e = 10, 12, 14, 16, 18 and 20 for - 0, 1, 2, 3, 4, 
and 5, respectively. The shape of the eigenfunctions was checked in 
order to avoid Calculations of higher modes. 
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CHAPTER IV 

NUMERICAL RESULTS AND DISCUSSION 

The results of stability analysis are usually presented as curves 
of Goertler number versus wave number for constant growth rate. Neutral 
stability curve (a = 0) separates the stable from unstable regions and 
increasing growth rate is akin to moving downstream with the vortices. 
Goertler number which represents viscous and curvature effects can be 
associated with a particular geometry and flow conditions. 

Effect of compressibility of the mean flow on Goertler instability 
is investigated for a range of Mach numbers from 0 to 5. These results 
are reported by El-Hady and Verma (1981a) and discussed in details in sec- 
tions 4.1 and 4.2. An attempt is made in section 4.1 to evaluate 
the effect of the mean density and mean viscosity on the stability charac- 
teristics of the boundary layer. The effect of using different thermal 
disturbance wall boundary conditions is examined in section 4.3. Effect 
of compressibility of the mean flow on the shape of the eigenfunctions 
is discussed in section 4.4, while its effects on the amplitude ratio 
of Goertler vortices is given in section 4.5. The effect of suction and 
cooling of the boundary layer on the development of these vortices is 
given in section 4.6 and 4.7, respectively. 

4.1 Effect of Compressibility on the Neutral Curve 

Figure 4 shows the neutral stability curves for different Mach num- 
bers. Instability of the boundary layer sets in at higher Goertler 
number as Mach number increases for all wavenumbers higher than 0.1. 

The critical value of G below which the flow is stable for any distur- 
bance wavenumber increases as Mach number increases. The stabilizing 



effect of compressibility is qualitatively in agreement with the results 
of Kobayashi and Kohama, but differs quantitatively due to the effect 
of the boundary layer growth that is included in the present analysis. 

It is worth noting that the incompressible limit, M*. = 0, of the neutral 
stability curve agrees with the results of Floryan and Saric and Ragab 
and Nayfeh. 

The rate of heat transfer between the fluid and the wall and the 
subsequent development of a thermal boundary layer in addition to the 
velocity boundary layer, play an important role in a compressible flow. 
The mean-density and the mean-viscosity variations inside the boundary 
layer have different roles in shaping the stability characteristics of 
the flow. Figure 5 shows the variation of the critical Goertler number 
with for a neutrally stable disturbance. The mean density is shown 
by curve 3 to locally destabilize the flow (curve 3 is calculated with 
constant mean viscosity). This effect can be predicted from Rayleigh 
criterion for invisqid flows with curved streamlines. His criterion 
for instability modified for a case with density variation is 
Py/p + (r Uco)y/(r Uoo) > 0 where r is the radius of curvature of 
the surface. This criterion shows that positive density gradient seems 
to destabilize the boundary layer. A compressible mean flow along an 
adiabatic flat plate has a positive density gradient as shown in Fig. 6. 
The figure shows the variation of the mean density and mean viscosity in 
the boundary layer for different Mach numbers. The effect of the mean 
viscosity is shown by curve 2 in Fig. 5 (calculated with constant mean 
density). The mean viscosity has a stabilizing effect compared to the 
mean-density variation. The exclusion of the disturbance viscosity, 
shown by curve 4, has a profound destabilizing effect on the 
boundary layer. It is worth noting that the conditions of 
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curve 4 are similar to the assumptions taken by Aihara in his calcula- 
tions. In spite of that, curve 4 indicates that critical Goertler 
number is increasing as Mach number increases which is opposite to 
that predicted by Aihara. It may be concluded that the stabilizing 
influence of the mean viscosity is more dominant than the opposite 
influence of the mean density, resulting in a more stable boundary layer 
as Mach number increases. 

4.2 Growth Rates of the Goertler Vortices 

The calculation of the initial point of instability or the minimum 
critical Goertler number as a function of mean flow or disturbance 
parameters is not that useful quantity regarding transition dependence. 
The growth of the vortices and not its initial point of instability is 
the decisive factor. It is one of the goals of this investigation to 
study the effect of Mach number on the growth rates of Goertler vortices 
and to present design charts for a range of Mach numbers. 

Figures 7-12 show contours of constant growth rates plotted in the 
G-3 plane for a range of Mach numbers from 0 to 5. The stable region 
below a = 0 curve is increasing as Mach number increases. The growth 
rate curves posses minimum which indicate a trend to higher wavenumbers 
at higher growth rates. They form a locus of the maximum growth rates 
of different wavenumber components. As Mach number increases, Figs. 7-12 
indicate that maximum growth rates occur at lower wavenumbers and higher 
Goertler numbers. The locus of the maximum growth rates shifts up and 
to the left of the chart as Mach number increases. 

Figure 7 (f^ = 0) and Fig. 10 (M^ = 3) show constant growth rate 
curves calculated by excluding all terms due to the growth of the 




boundary layer. These terms greatly influence the neutral stability 
curve specially at high Mach number. The effect of these terms dimin- 
ishes as the growth rate of the disturbance increases. 

For different Mach numbers. Fig. 13 shows the variation of the 
vortex growth rate with the Goertler number. The values plotted in 
Fig. 13 correspond to the minimum of the growth rate curves. The figure 
indicates that the vortex becomes more sensitive (its growth rate) to 
changes in Goertler number as Mach number increases. The influence of 
compressibility on reducing the growth rates is very small at high 
Goertler numbers. Fig . 13 also shows that compressibility has its 
maximum stabilizing influence when the vortex is weak. This conclusion 
is best illustrated by Fig. 14 where the parameter (G oc -G c )/G oc is used 
as ordinate to indicate the stabilizing effect (increase of critical 
Goertler number) as Mach number increases compared with G oc» the critical 
Goertler number at = 0 at the corresponding value of the growth rate a 


It should be reminded that the present definition of Goertler num- 
ber is based on the reference length L* defined in Eq.(l), whereas the 
original parameter identified by Goertler is based on the momentum thick- 
ness of the boundary layer. For the purpose of comparison with data re- 
ferenced to the displacement thickness 61 or the momentum thickness <$ 2 , 
the following relations are used: 

G n “ G L A n ' p n - B L“n 


3/2 e„ ■ 61 °n ’ °L A n 


where the subscript L denotes data based on the reference length L , 


and n denotes data based on the displacement thickness 61, or the momen- 
tum thickness 62 as reference length. The following table gives values 
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for the conversion factor A n for different Mach numbers, and freestream 
stagnation temperature of 310 K s where A-j = <^i /L , ^ ~ 



0 

1 

2 

3 

4 

5 

A 1 

1.720 

2.129 

3.365 

5.48 

8.504 

12.393 

A 2 

0.663 

0.656 

0.643 

0.641 

0.645 

0.646 


4.3 Effect of the Wall Boundary Condition of the Thermal Disturbance on 

The Stability Characteristics 

All the results presented so far are for an adiabatic wall. 

The wall boundary condition of the temperature disturbance requires a 
solution of the heat conduction equation very close to the surface. This 
leads to a boundary condition in the form of Eq. (21), that should be 
applied to a particular solid and a specified disturbance wavenumber. A 
detailed study of this boundary condition, with its two extreme cases 
of B ■* 0 (0y= 0, the wall is a perfect insulator), and B + (0 =0, 

the wall is a perfect conductor), leads to the conclusion that local 
stability near to the neutral region is dependent on the temperature 
fluctuation and the boundary condition imposed on it in a compressible 
boundary layer. At - 3, Fig. 15a shows a comparison between neutral 
stability curves calculated using the wall condition 0^(0) = 0 and 
$(0) =* 0. Apart from the differences in the value of the critical 
Goertler number, large discrepencies also appear for disturbances 
having low wavenumbers. The reason is that for low wavenumber distur- 
bances the thermal fluctuation can penetrate large distances into the 
solid wall, and hence, the wall temperature cannot remain at the wall 



temperature of the mean boundary layer. Figure 15b shows the effect of 
the thermal wall boundary condition on critical Goertler numbers of 
neutrally stable disturbances. The wall condition 0^(0) = 0 shows 
always a locally stabilizing effect in a compressible boundary layer 
compared to the condition 6(0) = 0, regarding the neutral stability 
curve. Using the full thermal boundary condition, Eq. 21 , the actual 
neutral stability characteristics will lie, probably, somewhere between 
the curves presented in Figs. 15a and 15b. It is of interest to notice 
that the effect of the thermal boundary condition vanishes with the 
growth of the disturbance. Figure 15c shows that for M^ = 3 which is 
typical for all the range of Mach numbers under investigation. 

4.4 Effect of Compressibility on Eigenfunctions 

Figures 16a-16d give a comparison of the shape of the eigenfunctions 
of u, v, w, and 0 , respectively at different Mach numbers for a distur- 
bance having a wavebumber 3 = 0.3 and a growth rate cr = 5. The cor- 
responding Goertler numbers are G = 13.4637, 12.9289, 12.0844, 11.5310, 
11.4786, and 11.7611 at = 1, 2, 3, 4, and 5, respectively. The 
values of u, v, w, and e are normalized with the maximum of the u com- 
ponent at the corresponding Moo. Figure 16 shows that the location of 

A A A A 

u , v • , w .. and 0 . move away from the wall as Mach number in- 
max min max min 

creases. The case of neutral stability (not shown) indicates a persis- 
tence of the disturbance outside the boundary layer that increases as 
Mach number increases. At H» = 3, Fig. 16a shows the shape of the 

A 

eigenfunction u when terms due to boundary layer growth are neglected 
for a neutrally stable as well as a growing disturbance of cr = 5. The 
effect of the boundary layer growth is to move the location of u max 
away from the wall leading to a slow decay of the disturbance outside 
the boundary layer. 
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In their study on the influence of suction or cooling on Goertler 
instability in compressible boundary layer's, El-Hady and Verma (1981b) 
showed that high suction or cooling rates bring the location of u |nax 
v m in> w max » ancl ®min nearer to the wa ^ ancl result in stabilizing the 
boundary layer. The stabilizing mechanism here is different from that 
due to compressibility. The effect of suction and cooling is further 
discussed in secs. 4.6 & 4.7. 


4.5 Effect of Compressibility on the Amplitude Ratio 

In this section, the growth of the measured by the amplitude 

ratio is considered since the growth rate alone is of little help 

in correlating transition. The logarithmic derivative of the amplitude 
of any disturbance quantity q in Eq. (15) is 

(In q) x = a (38) 

In terms of G instead of x and introducing a for q 
In (a/a Q ) = (4/3) /jjj (o/G) dG 

where a Q and G Q are the amplitude and Goertler number at the beginning 
of the unstable region, and the integration is carried out along a 
particular growth path. 

In all experiments of fncompressfbTe flows along concave surfaces, the 
growth path of the vortex was determined from the conservation of 
its dimensional wavelength A* that was observed in the flow direction. 

With lack of corresponding information for compressible flows and for 
comparison purposes, we assume the same criterion to hold. Since 
the dimensionless wavenumber 3 = 3*L* varies with x* for constant 
A*, the dimensionless wavelength 


A = 


U A 
00 


X ,1/2 




(40) 
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is used instead. In the G-8 plane, lines of A * const, are straight 
with slope 3/2. The locus of maximum growth rates closely coincide 
with these lines for the corresponding Mach number. We assume that a 0 
in Eqi39)is independent of A and G so that the dependence of a on A 
and G will be the same as that of a/a Q . For different Mach numbers, 
a series of amplitude ratio curves are calculated by integrating Eq. (39) 
for constant A up to G = 20. Figure 17 shows a reduction in the maximum 
value of the amplitude ratio as well as a shift of the most unstable 
wavelength to higher value as Mach number increases. Fig. 17 shows 
also an increase in the upper band of the disturbance wavelength that 
is always attenuated (cut off wavelength) as Mach number increases. 

It is a known fact that the stability theory cannot predict which 
disturbance wavelength will actually appear for a given surface con- 
figuration and flow conditions. For incompressible flows along con- 
cave surfaces, the wavelength of the distrubance appears to be deter- 
mined by the particular edge effects of the experimental apparatus 
(Tani and Sakagami , 1962) or by the oncoming disturbances (Bippes, 

1978). Based on incompressible experimental data (Tani and Sakagami 
1962; Tani, 1961), Floryan and Saric (1980) suggested that the wave- 
length selection mechanism (which decide the vortex growth path) may 
be based on the maximum growth rate of the disturbances. However, 
the selection process may be very easily affected by the properties of 
the experimental apparatus that determine the entry location to the 

locus of the maximum growth rate. For different Mach numbers, we 
integrate Eq. (39) along the locus of the maximum growth rate. The results 

are shown in Fig. 18 as function of the Goertler number. The integra- 
tion is carried up to G = 20 for comparison purposes. The figure 
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shows a stabilizing effect as Mach number increases in terms of lower 
values of ln(a/a 0 ). It is worth to notice that the values of ln(a/a Q ) 
at G = 20 for different Mach numbers following the path of maximum 
growth rates in Fig. 18 are almost indentical with the corresponding 
values at the most unstable wavelengths in Fig. 17 . 

4.6 Effect of Suction on the Stability Characteristics 

It, was shown by Goertler (1954 ) and Hammerlin (1955) that the curve 
of neutral stability depends insensibly upon the basic velocity profile 
in the laminar compressible boundary layer, when the changes in its 
momentum thickness is very small. However, Kahawita and Meroney (1977) 
showed that the inclusion of the normal velocity terms change the 
location of the neutral curve drastically. In section 4.2, the same 
conclusion was reached for the compressible boundary layer. Terms due 
to boundary layer growth have large effect near the neutral stability 
region specially at high Mach numbers. Therefore, it is expected that 
variation of the velocity profile due to suction (which does not affect 
the momentum thickness but changes the normal velocity term) may change 
the critical Goertler number. 

The effect of suction of the laminar boundary layer on Goertler 
instability is examined at M*, - 0.8 and 3. The self similar suction 
parameter y defined in Eq. (5) is used for this purpose. Figures 19 and 
20 show the change in the location of the neutral stability curves for 
different values of the suction parameter at 0.8 and 3 respective- 

ly. It is observed that critical Goertler number first decreases with 
increasing suction (destabilizes the boundary layer). The boundary layer 
with a suction given by y & -0.45 and -1.3 is the most unstable at 
M ot = 0.8 and 3, 


respectively, as far as the critical Goertler number is concerned. By 
increasing the suction parameter above these levels for each Mach num- 
ber, the critical Goertler number increases. Fig. 21 shows how the 
critical Goertler number changes with y for different Mach numbers. The 
level of suction that is required for a noticeable increase of the 
critical Goertler number increases as Mach number increases. This 
suction level is of higher order of magnitude than the suction required 
to stabilize Tollmien-Schl ichting waves (see Saric and Nayfeh, 1977 for 
a 0). This leads to the conclusion that with these levels of suction, 
Tollmien-Schl ichting waves are practically eliminated, and Goertler 
vortices may dominate the flow. 

At = 3, Figs. 10, 22, and 23 show waves of constant growth rates 

for y - 0, -1.2 and -1.6. By increasing suction, curves of small growth 

rates are compressed and moved to higher Goertler numbers. Curves of 
high growth rates are slightly affected by suction especially at high 
Mach numbers. 

To see the overall effect of suction, the growth of the vortices 
should be taken into account as they develop downstream. For M = 3 
Fig. 24 shows the amplitude ratio of the vortices as a function of 
Goertler number. Integration of the growth rates is performed using 
Eq. (39) to G = 20 for comparison purposes. It is clear that despite 
indications of the destabilizing effect of small suction for all Mach 
numbers regarding the critical Goertler number (see Fig. 20), the 

overall effect of suction is to stabilize the boundary layer as shown 
by Fig. 24 for = 3 and by Floryan and Saric for = 0. Although 

a suction corresponding to y = -1.2 which is near to the critical value 
of the suction parameter at M ro = 3, locally destabilizes the boundary 
layer in terms of the critical Goertler number, it 
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reduces the amplitude ratio by about twenty- two percent thus sta- 
bilizing the boundary layer. 

At different Mach numbers. Fig. 25 displays the amplitude ratio 
of the vortices at G = 20 over a range of the suction parameter y. The 
amplitude ratio of the vortices reduces as suction increases. With 
increasing Mach number, it becomes more difficult to stabilize the 
boundary layer (reducing the amplitude ratio of the vortices) unless 
very high suction is used. Figure 25 shows that at = 5 the ampli- 
tude ratio is hardly influenced by suction in the range of 0<y<2. 

At Mo, " 3, Fig. 26a-26d show a comparison of the shape of the 
eigenfunctions u, v, w, and 0, respectively, for different values of 
the suction parameter y( normalized with the maximum of u component). 

These eigenfunctions are for a neutrally stable disturbance having a 
wave number 3 = 0.3. The corresponding Goertler numbers are G = 

1.1155, 1.7549, 3.5193 at y = -0.6, -1.6, -2.0, respectively. 

ys /\ ^ 

Figure 26 shows that the location of u , v . , w , and 0 . 

max min max min 

move towards the wall as suction increases. 

The normal velocity component, if directed away from the wall (the case 
of no suction), tends to destabilize the flow by encouraging penetration 
into the free stream where viscous dissipation is small. By increasing 
suction the thickness of the boundary layer decreases, and the distur- 
bance is confined to a region closer to the wall where dissipative 
action Is strong, thereby Increasing the stability of the flow. 

4.7 Effect of Cooling on the Stability Characteristics 

At Moo = 3, Fig. 27 shows how the location of the neutral curve varies 
with different cooling rates. The cooling parameter 0 W / G ad is used for 
this purpose, where ^,/ Q ad ~ 1 represents no cooling, and V 0 ad < 1 re ‘ 
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presents cooling. Small cooling has a destabilizing influence 4 , however 
as cooling increases, the neutral curve moves upward indicating a sta- 
bilizing influence as far as the critical Goertler number is concerned. 
For wavenumbers approximately greater than 0.4, Fig. 27 shows that cool- 
ing always destabilizes the boundary layer (Goertler number decreases). 

It should be noted here that in section 4.3 we arrived at the con- 
clusion that neutral stability curve is dependent on the wall boundary 
condition of the thermal disturbance. Therefore, one may think that 
comparing neutral curves in case of cooling, where the wall is consi- 
dered completely conducting and the wall boundary condition of the ther- 
mal disturbance is 0 = 0, with that in case of no cooling, where the wall 
is considered ideally insulated and the wall boundary condition of the 
thermal disturbance is 0y = 0, is somewhat misleading. The neutral sta- 
bility curve is calculated for a no cooling case with the thermal boun- 
dary condition 0 - 0 at y = 0, and is shown in Fig. 27 for comparison. 
Previous conclusions regarding effect of cooling on the critical Goertler 
number are still valid. 

At = 3, Fig. 28-30 show contours of constant growth rates for the 
cooling parameter 0^/O a( j = 0.75, 0.25, and 0. 15, respectively. These 
figures display the influence of cooling on the vortices at different 
stages of their growth. Condi usions are summarized in Fig. 31, where the 
parameter (G nc - G c )/G nq is used as ordinate to indicate the stabilizing 
or destabilizing effect (increase or decrease in critical Goertler num- 
ber) as cooling increases compared with G n c> the critical Goertler num- 
ber at 0w/O a d * 1 for the corresponding values of the growth rate a. 

Fig. 31 shows that at Moo B 3, weak vortices (small growth rates) are more 
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influenced by cooling. Cooling parameters 1 > 0 ^/ 03 ^>0. 75 seems to de- 
stabilize the boundary layer. As cooling increases (0 w /0 ac | < O.75) , Fig. 

31 shows a stabilizing influence of cooling on weak vortices; however, 
cooling has almost no influence on strong vortices (high growth rates). 

The above conclusions represent a local effect of cooling on the sta- 
bility limit as well as growth rates of the vortices. It is important 
to estimate the effect of cooling on the total growth of the vortices 
by calculating the amplitude ratios. Again, Eq. (39) is used to in- 
tegrate the growth rates along the locus of maximum growth to G = 20. 

At M*, = 3, Fig. 32 shows a comparison of the amplitude ratios for 
different cooling parameters. The figure shows that stabilizing the 
boundary layer with respect to Goertler vortices can hardly be achieved by 
using very high cooling. At G = 20, the amplitude ratio of the vor- 
tices reduces by about 2 % at - 0.25 and by 5.6 % at ©w/Oad = 

0.15. The case of 0 w /0ad = 0.75 shows a destabilizing effect at this 
Mach number. 

Fig. 33 shows the effect of cooling on the amplitude ratio of 
Goertler vortices for different Mach numbers. The figure clearly dem- 
onstrates that moderate cooling has almost no influence on the stabil- 
ity. However, high cooling influences low Mach numbers more than high 
Mach numbers. At O w /O ac | - 0.25, the amplitude ratio is reduced by 13 . 4% at 
Ko - 0.8, by 2.8% at Mo = 3, and increased by 2.6% at Mco = 5. 

It is clear from Figs. 32 and 33 that the effect of moderate cool- 
ing is relatively small considering the strong temperature dependence 
of the viscosity and density and the corresponding large variation of 
both in the boundary layer. At M* a 3, Fig. 34 shows the distribu- 
tion of viscosity and density of the mean boundary layer flow for differ- 



ent cooling parameters. In contrast to the insulated wall, where the 
density gradient is positive and the viscosity is large at the wall, 
the densi-ty profile in the cooling case has not only positive but neg- 
ative gradients. The viscosity profile has also two regions where 
P > 1 and P < 1. The negative density gradient may have a stabilizing 
effect according to the Rayleigh inviscid criterion, but the varia- 
tion of both density and viscosity in the cooling case are complica- 
ted enough to compare with the case of no cooling. 

The small effect of moderate cooling shown in Figs. 32 and 33 
is in marked contrast to the results for the Tol lmien-Schl ichting (TS) 
type instability which is extremely sensitive to wall temperature 
(see Fig. 7 in Mack, 1975). This is ultimately due to the presence 
of an inner critical layer in TS disturbances through which the eigen- 
function varies rapidly. Therefore, stability characteristics depend 
critically on the local properties of the mean-flow velocity, viscosity, 
and density distributions which are very sensitive to the wall tempera- 
ture. For Goertler instability, however, there is no inner critical 
layer and the centrifugal force is the controlling factor. This insta- 
bility depends only on the overall properties of the mean flow, such 
as the average velocity gradient, which are much less influenced by 
the wall temperature. 

At Mo * 3, Figs. 35a-35d give a comparison of the shape of the 
eigenfunctions of u, v, w, and §, respectively for different values of 
the parameter 9 w /O ac j* These eigenfunctions are for a neutrally stable 
disturbance with wavenumber 3 = 0.3. The corresponding Goertler num- 
bers are G -= 1.463, 1.301, 1.205, and 1.213 for 0 w /0 a d = 1, 0.75, 0.50 
and 0.25, respectively. The values of u, v, w, and © are normalized 
with the maximum of u - component at the corresponding cooling 
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parameter. The figure shows that the location of u max , w max , and 
§ ml * n move towards the wall as cooling increases. The velocity compo- 
nent v shows a persi stance outside the boundary layer. 

Comparing the shape of eigenfunctions due to cooling in Fig. 35, 
with that due to suction in Fig. 25, shows that at Moo = 3, both high 

suction and high cooling stabilize the boundary layer by con- 
fining the disturbance to a highly dissipative region nearer to 
the wall. Suction may be more effective than cooling as it brings 
the disturbances more close to the wall. 
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CHAPTER V 
CONCLUSIONS 

A leading order approximation of a compressible stability theory 
of boundary- layer flow along a concave surface is presented and solved 
numerically. The compressible boundary layer equations along a flat 
plate are used to represent the mean flow. The results of the stability 
analysis are summarized as follows: 

1. Instability of the boundary layer with respect to Goertler 
vortices sets in at higher Goertler number as Mach number 
increases. 

2. The local stability near to the neutral region is dependent on 
the wall boundary condition imposed on the temperature fluctu- 
ation in a compressible boundary layer. 

3. At high Mach number, the growth of the vortices is sensitive 
to small changes in Goertler number. 

4. Compressibility has its maximum influence on the vortex when 
it is weak. 

5. Terms due to boundary layer growth have large local effect 

near the neutral stability region specially at high Mach numbers. 

6. Compressibility reduces the maximum amplitude ratio by about 
20% as Mach number increases from 0 to 5. 

7. With Increasing Mach number, the most unstable and cut off 
wavelengths shift to higher values. 

8. Suction may have a local destabilizing effect on the boundary 
layer as far as the critical Goertler number is concerned, if 
its level is below a critical value. This critical level of 
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suction, where local effect becomes stabilizing increases 
with Mach numbers. 

9. In contrast with local suction effects, suction of the boun- 
dary layer has always a stabilizing influence on the total growth 
of Goertler vortices. The amplitude ratios of the vortices 
reduce with any level of suction used. 

10. Suction of the boundary layer is more effective in stabilizing 
the flow at low Mach numbers than at high Mach numbers. It 
becomes increasingly difficult as Mach number increases to 
reduce the amplitude ratios of the vortices unless very high 
levels of suction are used. 

11. Wall cooling, like suction, has a local destabilizing effect 
on the boundary layer as far as the critical Goertler number 
is concerned. 

12. At low Mach numbers, small cooling applied to the wall 
(Ow/Oad — 0.75) has no influence on the amplitude ratios 

of Goertler vortices. A noticeable reduction in the ampli- 
tude ratio starts with moderate cooling (0w/Oad £0.5). 

13. At high Mach numbers, it seems almost impossible to stabilize 
the boundary layer using practical rates of wall cooling. 

The amplitude ratios of the vortices increase and the boun- 
dary layer becomes more unstable with small or moderate cooling. 

14. Goertler instability is more difficult to influence and control 
by suction or cooling than Tollmien- Schlitching instability. 

The reason is that Goertler instability which is referred to 

as centrifugal type instability, depends on the velocity differ- 
ence between the inner and outer region of the boundary layer 
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and not on the shape of the boundary layer profile as in the 
case of Tollmien-Schl itching instability. 
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Figure 1 Body oriented coordinate system. 



Figure 2 Coordinate system based on streamlines and 
potential lines. 
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Figure 5 Effect of mean density, mean viscosity, and 
disturbance viscosity on minimum critical Goertler 
number, in a compressible boundary layer along an 
adiabatic flat plate. 






Figure 7 Contours of constant growth rates 
at Mach number = 0. Terms due to boundary 

growth Included , terms due to boundary 

layer growth excluded aary 



Figure 8 Contours of constant growth rates 
at Mach number - 1 . 
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Figure 10 Contours of constant growth rates 
at Mach number = 3. Terms due to boundary 

layer growth included , terms due to boundary 

layer growth excluded 
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Figure 14 Effect of compressibility on local 
stability. 
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Figure 15 Effect of the wall 
boundary condition of the thermal 
disturbance on the stability 
character!* stfcs. 










Figure 16 Effect of compressibility on che shape of the eigenfunctions for a 
disturbance having a wavenumber 0.3 and a growth rate 5. Terms due to boundary- 
layer growth included , terms due to boundary-! ayer growth excluded 
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Figure 18 Effect of compressibility on 
the amplitude ratio calculated along the 
locus of maximum growth rates. 











Figure 22 Contours of constant growth rates at and 


suction parameter y=-1.2. 





Figure 23 Contours of constant growth rates at M ro =3 and 


suction parameter y=-1.6. 
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Figure 25 Effect of suction on the amplitude ratio at 
G=20 for different Mach numbers. 
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Figure 26 Effect of suction at on the shape of the eigen- 
functions for a disturbance having wavenumber 0.3 and 
zero growth rate. 
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Figure 27 Neutral stability curves for different cooling rates 

at M =3. 
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Figure 28 Contours of constant growth rates at M =3 and 
0w/0ad=O-75. 







Figure 30 Contours of constant growth rates at M =3 and 

oo 

®w/®ad = ^* ^ 
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Figure 34 Variation of mean density and viscosity for different 


cooling rates. 
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APPENDIX A 

This appendix contains the dimensional field equations for com- 
pressible flow (see Hughes and Gaylord, 1964) written in the specified 
system of coordinates x, y, z. The velocities are given by u, v, w, 
and the metric coefficients are h, h, 1 in the x, y, z directions re- 
spectively. The Prandtl number r and specific heat Cp are taken con- 
stant in the stability analysis. The ratio c of the second to the first 
mean flow viscosity coefficients is defined as c = 2(e - l)/3, e is 
taken equal to 0.8. 


x-momentum 


pjX + hV wu z ' ^ vh x ' uh y J ] = ‘ K 

+ F (llV,v) X + ^[ {Zllh( K + 7 h y )} x 
+ tp h [<F>x<Jy {yh 2 ( vK )} z] 
[^x + (F)y] h y-^[F v y + 7 h x] h : 


+ JL 
h 2 


[F v x + F v y + wv z + 7 (vh x- ljh y ) ] = ‘F p y 
+ f(liV • v) y + ;r[{uh [(^x + ( F } y] } x 

+ ™ (F v y + ^ h x } } y + ^ ( F w y + v z } } z] 

[<F>x + <F>y] 


+ \ 
n 




y-momentum 



z-momentum 


Energy 


Continuity 


P[F w x + F w y + WW 2 j = * p z + c(pV * v) z 
+ ^[{uh(u 2 +iw x )} x + tph(i-w y + v z )} y 

+ {2Mh 2 w z ) z J 

pc p^F °x + F Q y + w 0 z^ ” F p x + F p y + wp z 

+ < c + 2 > p [(f u x + 7V Z + { F v y + JV 2+w z 

+ p [^ w y + v z ) 2 +(u z + FT w x ) 2 +{( F ) x +( F ) y }2 
+ 2c((i vJV ( V/Jr h x> ,( i U x 
+ 7 h y )w z +( F v y + 7 h x )w z}] 

+ ^[ (u 0 x ) x + (w 0 y ) y + h 2 (|J 0 z ) z] 

(hpu) x + (hpv) y + h 2 (pw) z = 0 

where 

V • v --V[(hu) x + (hv) + h 2 w z ] 



APPENDIX B 


The two-dimensional compressible boundary layer equations for 
zero pressure gradient together with the equation of state for a perfect 
gas are reduced to the following set of ordinary differential equations. 

(uii y ) y + gu y - YU y (bi) 

( r Vy + 9G y + 2)j( U y )2 * ^ 0 y = 0 (B2) 

g y - i pU = 0 (B3) 

U = 0, g = 0, and 0^ = 0 or 

0 = 0 W at y = 0 (B4) 

U -> 1, 0 + 0 as y + y e (B5) 

by using the transformation 

y = (U*/v* (p*/p*) dy* (B6) 

as well as a stream function to satisfy the continuity equations. Here 
0 = (i* - ie)/(ioe ” O’ w ^ ere io are ^e fluid enthalpy and stag- 
nation enthalpy respectively, and e denotes conditions at the edge of 
the boundary layer. Equations B1-B5 are integrated numerically with 
the thermodynamic and transport properties of the perfect gas compu- 
ted at each integration step. The variation of perfect gas properties 
with temperature are taken from (Hilsenrath, Beckett, et a!., 1955). 
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APPENDIX C 


This appendix contains the nonzero elements of eigenvector matrix 
D for the case when X 2 anc * ^6 are re P eate d* 


d 3J 

d 4J 

Dyj 

d 8J 
J = 2,6 

°3J 

°7J 

d 8J 
J = 3,7 

D 1J 

°2J 

d 3J 

D 4J 

°7J 

d 8J 
J = 4,8 

d 3J 

d 40 

d 5J 

°60 


-B/Aj 

8[o + B 2 + Aj(V e - Aj)]/Aj 
1.0 
X J 


-B/Aj 

1.0 


B(V e - 2 Aj) ( 1 - Aj/ 8 2 )/(2G 2 + V xe )Aj 

X J D 1J 

1.0 

(V e - 2Aj)/B 

1.0/Aj - Aj/8 - ctO^j/B 

1 + x J D 7J 


W + W - $1 - a^ ta^t 

-Xj[o + $ 2 + A j(V e - Aj)]D 3J /B 2 

1.0 


J 

-AjD 3J /B + (a + XjV e )/B 


d 8J = x J d 7J 
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APPENDIX D 

This appendix contains the nonzero elements of the eigenvectors 
matrix D for the special case a = &V e whereXi, X2, ^ 3 » and M are 
repeated. 

J = 1,5 

° 7 J = -\)/B 
° 8 j ~ ’3 
J = 2,6 

DlO " 1 
° 2 J = \j 
D 3 J = 1 
D 5 J = 2 

d 6J = 2x J 
d 7 J = -\}/3 
D 8 j - -6 
J = 3,7 

° 3 J " 1 

D., = 2 
4 J 

D7J - -(1 + Xj)/8 

D 8J = + 2 ^ /(3 

J = 4,8 

d 3J a 1 

d 4 J = *j( 2 \j + 3 )/ 3 2 
D 5J * 4 /G 2 
D 6 J = 4 Xj/G 2 
D7J = -(1 + Xj)/8 


d 8j = -0 + x j) 2 /b 



APPENDIX E 


This appendix contains the nonzero elements of the eigenvectors 
matrix for the Adjoint Problem (34). 

0 = 1,3 

DlJ =' i v[f-(aV e - 2G 2 - V xe ) - 0 ( a + B 2 ) * (a + B 2 )D ? , ] 

0 o 3 

• %*"*' - 261 - v - 

&3J ' (o + S 2 )/$ 


d 4J = B/Aj 

d 5J = - |^t(are +g 2 )D6j + (c + 2) c g + ^ + ^(G 2 + V xe - V e , 
(c + 2)V e (ro + B 2 )] 




- - BV - + (q + 2 ) 


e x 


r e v e (v e e - ^)] 


D7J = -v e 

J = 3,7 

D u = ^t( a +s2 HxJ - 1-0) + (°\> + 2G 2 + v xe )o 4J 
D 2 j - 1.0 
d 3 j = 


'4J 


Xj “4J 

IVfi+iAlL. 

(2G 2 + Vxe) 


°5J = -^jCC r e cr +e 2 )D 6J + [(c + 2)V e (r e a +B 2 ) + G 2 + V xe + 
(q + 2)\]a]D4j] 


d 6J = ID 2 — V^P - (c 1 2)Vp (r P V 0 + X? - r „0 -B 2 np A .i 
(r e v e + Xj 2 . r e a -b 2 ) 
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nnnnooononnonnnn 


APPENDIX F 


FORTRAN PROGRAM FOR STABILITY' ANALYSIS. 

. THIS PROGRAM CALLS SUBROUTE JOCK (LINE A 430) WHICH 
IS A VARIABLE STEP SIZE INTEGRATOR WRITTEN BY SCOTT AND WATTS (1977). 


PROGRAM STATIC (INPUT , OUTPUT* TAPES* INPUT , TA PF 6-CUTPUT, TA PF 10 ) 

REAL K ( 8* 8 ),KK (8, 8),KNT,KNTB,KTB,MACH,MU,MUP 

DIMENSION Z ( 3, 101 ) » A(4,P), AL PHA ( 4 ) , 8(4,8)# BETA ( 4 ) , K ARR AY( 7 ) , 

1 C ( 8 , 8 ) , CCC ( 8, 8 ) , FF ( 8 , 8 ) , DD(8,8), XL C 8 ) 

DIMENSION WOP K ( 10000 ) , IW0RK(150) 

DIMENSION P V ( 8 ) , IK(30), C0(8) 

COMMON / AAA/ KNT , B E T, V EL F , GG, M AC H, CC , P X , P Y 
COMMON ABBB / XS AVE , KL , I NDE X 

COMMON /CCC/ V(101 ), VP ( 101 ), VP P( 101), DUX (101 ),DUPX(101), DVXC101 ),D 
ITXC 101 ),DTPX( 101) ,DMUX (101 ),DALFX( 101 ) 

COMMON 7000 / Y ( 101 ) ,U ( 101 ) ,UP ( 101 ) ,UPP ( 101 ) , T ( 101 ) , T P ( 101 ) , TPP ( 101 
1 ),PRANOL(101),MU(101),MUP(101),ALFA(101),ALFAP(101) 

COMMON /FFF/ IE, MOD 
COMMON /EEE/ ACC, FACT 

NAMELIST /GORTLR/ GG, KNT , BFT, AC C , MOO, I TR , I NDE X ,N IT, I L M, F ACT 
******************************************************** 

MOD INDICATES DIFFERENT MODELS 

MOD* 1 WITHOUT TERMS DUE TO BOUNDARY LAYER GROWTH 

MOD-2 WITH TERMS DUE TO BOUNDARY LAYER GROWTH 

ITR INDICATES ITERATION APPLIED TO 

ITP-1 GG KNT -CONST ANT 

IT»-2 KNT GG-CONSTANT 

INDEX INDICATES TYPE OF PROBLEM 

INOEX-l HOMOGENEOUS PROBLEM 

INOEX-2 ADJOINT HOM. PROBLEM 

NOTE- INDEX-1,2 ARE VALIO FOR MOD-2 ONLY 

NIT-MAX, NO, OF ITERATION 

IN IS NO OF ITFR ATIONS 

I M • 0 CONVERGENCE 

IE-NO OF POINTS IN THE Y DIRECTION 


A 1 
A 2 
A 3 
A 4 
A 5 
A 6 
A 7 
A 8 
A 9 
A 10 
A 11 
A 12 
A 13 
A 14 
A 15 
A 16 
A 17 
A 18 
A 19 
A 20 
A 21 
A 22 
A 23 
A 24 
A 25 
A 26 
A 27 
A 28 
A 29 
A 30 
A 31 


CD 

UD 


on on ooonoooo 000000000000000000 


vo 

o 


IM-1 NO CONVERGENCE 

It IS COUNTFR FOP THE NO OF POINTS ON PLOT 

ILf* IS MAX NO OF POINTS ON THE PLOT 

GG-GOERTLEP NUMBER 

KNT* WAVE NUMBER 

MACH* MACH NUMBER 

MU « VISCOSITY 

XL ( I )• EIGENVALUE ARRAY 

C ( 8, 9 ) “CONSTANT COEFFICIENT MATRIX 

K(8,8)-EIGENVECT0R MATRIX 

BET -GROWTH PATE 

VELF -NORMAL VELOCITY COMPONENT OUTSIDE THE B.L. 

CC -CONSTANT-O. A/3, 

ACC -ACCURACY 

FACT -A PARAMETER WHICH DETERMINES THE STEP SIZE 

IN NEWTON-R A PHSGN ITERATIVE METHOD. SEE SUBROUTINE 
ITRN. 

SUFFIX P DENOTES DERIVATIVES WITH RESPECT TO Y 
COORDINATE. 

SUFFIX X DENOTES DERIVATIVE WITH RESPECT TO X 
COORDINATE. 

THUS UP AND UPP ARE FIRST AND SECOND DERIVATIVES 
WITH RESPECT TO Y AND DVX AND OTX ARE THE X 
DERIVATIVE OF THE NORMAL COMPONENT OF VELOCITY 
AND TEMPERATURE 
PR ANDL ■ PR ANDTL NUMBER 

ALFA -DERIVATIVE OF MF AN VISCOSITY WITH MEAN 
TEMPERATURE. 

****************************************************** 

READ (10) ETA, MACH, IE 
READ (10) (Y(I ),I»1,IE) 

READ (10) (U(I ),I»1,IE) 

READ (10) (UP(I),I-1,IE) 

READ (10) (UPP ( I ) , I «1 , I E ) 

READ (10) (TCI), I -1, IE) 


A 32 
A 33 
A 34 
A 35 
A 36 
A 37 
A 38 
A 39 
A 40 
A 41 
A 4? 
A 43 
A 44 
A 45 
A 46 
A 47 
A 48 
A 49 
A 50 
A 51 
A 52 
A 53 
A 54 
A 55 
A 56 
A 57 
A 58 
A 59 
A 60 
A 61 
A 62 
A 63 
A 64 
A 65 
A 66 
A 67 


non 


READ (10) (TP ( I )» I«l> IE ) 

READ (10) (TPPCI)#I-1>IE) 

READ (10) ( PR ANDL ( I )> I»1 * IE ) 

READ (10) (MU(I)>I*1*IF) 

READ (10) ( MUP ( I ) > I *1 > IE ) 

READ (10) (AIFA(I),I-1,IE) 

READ (10) (ALFAP(I),I*1,IE) 

READ (10) ( DUX € I ) * 1*1 » I E ) 

READ (10) (DUPXU ),I*1,IF ) 

READ (10) ( DTX ( I ) * I * 1 > I E ) 

READ (10) ( DTP X ( I ) » I ■ 1 » I F ) 

READ (10) ( DMUX ( I ) * I*l> I F ) 

READ (10) (DALFX(I),I«1,IF) 

READ (10) (V(I)#I«1#IE) 

READ (10) (VP(I)f I»lfIF) 

READ (10) ( VPP ( I ) * I *1 * I E ) 

READ (10) ( DVX (I)»I a l»IE) 

C********* ******************************************* 


INPUT PARAMETERS 

£************************************************** 
PX-0.0 
PY-0.0 
CC—O. A/3.0 
CCl-COl. 

CC2-CC+2. 

PR1-PRANDL (1) 

VEL F*V ( 1 ) 

IL-0 

READ (5,G0RTLR) 


A 68 
A 69 
A 70 
A 71 
A 72 
A 73 
A 7 A 
A 75 
A 76 
A 77 
A 78 
A 79 
A 80 
A 81 
A 82 
A 83 
A 8 A 
A 85 
A 86 
A 87 
A 88 
A 89 
A 90 
A 91 
A 92 
A 93 
A 9A 
A 95 
A 96 
A 97 
A 98 


1 


IF ( EOF ( 5 ) ) A 5 » 1 
CONTINUF 

WRITE ( GORTLR ) 


A 99 
A 100 
A 101 



2 

CONTINUE 

A 102 


IF (AB$(BET)#EO.ABS(KNT*VELF) ) STOP 

A 103 


IN-0 

A 10* 


WRITE ( 6# 56 ) 

A 105 

3 

CONTINUE 

A 106 


IM-0 

A 107 

c********* ************************************** 

A 108 

c 

PARALLEL AND NON PARALLEL BRANCHING 

A 109 

c********************* ************************** 

A 110 


IF (MOD.NE.l) GO TO 5 

A 111 


DO * I *1# I E 

A 112 


V(I)»0.0 

A 113 


VP(I)-0.0 

A 11* 


VPP(I)-0.0 

A 115 


Dvxm-o.o 

A 116 


DUXUJ-O.O 

A 117 


DUPX (I ) *0* 0 

A 118 


DTXCD-0.0 

A 119 


DTPX ( I ) -0.0 

A 120 


DMUX ( I ) -0 • 0 

A 121 


DALFXII )-0.0 

A 122 

* 

CONTINUE 

A 123 


VEL F-V ( 1 ) 

A 12* 

5 

CONTINUE 

A 125 

c ******** ************************************** 

A 126 

c 

CONSTANT COEFFICIENT MATRIX 

A 127 

q*********** ************************************ 

A 128 


DO 6 1-1,8 

A 129 


DO 6 J-1,8 

A 130 

fc 

C(I# JJ-0.0 

A 131 


C ( 1, 2 ) -1 • 0 

A 132 


C(2,1)-KNT**2+BET 

A 133 


C(2*2)-VELF 

A 13* 


C(3#l)«- BET 

A 135 


C(3,5)-BET 

A 136 


C ( 3, 6) -VELF 

A 137 


C (3*7)— KNT 

A 138 



C(4, 1)— 2«0*GG**2*BET*VELF-DVX<1) A 139 

C(4,2)*-BET A 140 

C( 4,3) ■-KNT**2-PET A 141 

C(4,5)«GG**2+DVX(1)-VFLF*BFT+CC2*VFLF* (PR1*BFT*KNT**2) A 142 

C(4,6)*CC2*(BET+PR1*VELF**2)-VFLF**2 A 143 

C(4,7)«KNT*VELF A 144 

C ( 4, 6 ) *-KNT A 145 

C ( 5,6) *1 • 0 A 146 

C<6»5)«KNT**2+BET*PR1 A 147 

C ( 6, 6) «PR 1*VEL F A 148 

C ( 7, 8 ) »1 • 0 A 149 

C18,4) — KNT A 150 

C(8,5)»CC1*KNT*BFT A 151 

C(8,6)«CC1*KNT*VELF A 152 

C<8,7)»KNT**2+BET A 153 

C ( 8# 8 ) «VELF A 154 

IF (WOD.EO.l.AND.BET.EQ.OtO) GO TO 20 A 155 

KNTB-(VELF**2+4.0*<BET+KNT**2) )**0.5 A 156 

KTB*f(PRANDLCl)*VELF)+*2+4 # 0*(BET*PRANDL(l)*KNT**2))**«5 A 15 7 

C*********************************************************** A 150 

C EIGENVALUES A 159 

C *********** ****** ************ ******** ************ ********** a 160 

XL<1>— KNT A 161 

XL ( 2 ) *— 0« 5*(— VELF+KNTB) A 162 

XL<3)«-0.5*(-VELF+KNTB) A 163 

XLC4>— 0,5*CPRANDL Cl )* (-VELF ) +KTB ) A 164 

XL ( 5 ) *KNT A 165 

XL(6)»0»5*{VELF+KNTB) A 166 

Xl<7)»0.5*<VELF+KNTB) A 167 

XL (8)«0.5MPRANDL (1 )*VELF+KTB) A 168 

IF C INDF X • EO • 2 ) GO TO 13 A 169 

PS«2,0*GG**2tDVX(l) A 170 

PP«<GG*KNT)**2+DVX(1)*KNT**2 A 171 

C A 172 

DO 12 J*l, 8 A 173 

XLL-XL( J)4XL( J) A 174 

GO TO <7,8,9,10,7,8,9, 10), J A 175 


VO 

GO 



'X) 

■P' 


C*********** ********************************************** A 176 

C EIGENVECTOR MATRIX ELEMENTS A 177 

C****************** ******************* ************ ******** a 178 

C A 179 

7 K(l,J)-0.0 A 180 

K(2»J)«0.0 A 181 

K(3*J)»-KNT/XL(J) A 182 

K (** J )«KNT*( (BET+KNT**?) +XL( J )+( VELF-Xl ( J) > ) /XU A 183 

K(5*J)»0.0 A 184 

K(6*J>-0.0 A 185 

K(7,J)«1.0 A 186 

K(8#J)-XL(J) A 187 

GO TO 11 A 188 

8 K ( 1 # J ) «0, 0 A 189 

K<2>J)-0.0 A 190 

K(3#J)— KNT/XL( J) A 191 

K(4»J)«0.0 A 192 

K<5> J)«0.0 A 193 

K(6*J)«0.0 A 194 

M7, J)«1.0 A 195 

K ( 8# J ) «XL ( J ) A 196 

GO TO 11 A 197 

9 Xl«(KNT*(VELF-2.0*XL(J))*(1.0-(XL(J)/KNT)*+?))/(PS*XL(J)) A 198 

K(1»J)>X1 A 199 

K<2, J)-XL< J)*X1 A 200 

K(3> J ) -1 . 0 A 201 

M4,J)«(VELF-?.0*XL( J ) )/KNT A 202 

K ( 5 * J ) *0 • 0 A 203 

R(6»J)*0.0 A 204 

M7# J)-1.0/Xl< J)-Xl(J) /K NT-BET* XI /KNT A 205 

K(P,J)«1.0+XL(J)*K(7,J) A 206 

GO TO 11 A 207 

10 X4«BFT+KNT**2+XL(J)*<VELF-XL<J)) A 208 

K(1#J)*0#0 A 209 

K(2#J)-0.0 A 210 

K(3#J)»<XL(J)*(BET*XL(J)*VELF)-<PP/X4)) / ( X L ( J ) * * 2-KNT** 2 ) A 211 



K(4> J)—X4*XL(J)*M3> J)MKNT+*?) + (BET*XL(J)*VELF)*UCCU.0m4/KNT A 212 
1**2) A 213 

K ( 5# J ) *1 • 0 A 214 

K(6>J)»XIU) A 215 

K<7> J)«-XL( J )*K<3, J ) /KNTHBET + XL{J )*VELF)/KNT A 216 

K(8*J) a XL(J)*K(7*J) A 217 

11 CONTINUE A 218 

12 CONTINUE A 219 

GO TO 29 A 220 

C A 221 

13 CONTINUE A 222 

C******** ********************************* **************** A 223 

C AO JOINT PROBLEM E IGFN-VEC TORS A 224 

C***************************** ******** ****************** A 225 

BKN«BET+KNT**2 A 226 

BPK«BET*PR1+KNT**2 A 227 

DO 19 J*l» 8 A 228 

GO TO U4*15#16>17#14#15#16#17), J A 229 

C A 230 

14 CONTINUE A 231 

K(2f J)»(KNT*(BET*VELF-2.0*GG**2-DVXC1) ) /XL C J )-BET**2/KNT) / ( XL ( J >* V A 232 

1ELF-BET) A 233 

K(l#J)«-( (KNT/XL ( J ) )+(BET*VElE-2.0*GG**2-DVX(l ) )-BET*BKN/KNT+BKN*K A 234 
1(2*J))/XLU> A 235 

K C 3# J ) "BKN/KNT A 236 

K(4>J)«KNT/XL(J) A 237 

K(6# J)»<XL( J>*t 2.0*VFLE*BET/KNT-GG**?/KNT-DVX(1)/KNT)-KNT*VELF**2- A 238 
1 BET**2/KNT+CC2* PR 1*VELF* ( VELF*K NT-BET* XL < J )/KNT) ) / { PP1* { BET-VELF+X A 239 
?L(J))> A 240 

K<5> J)— (BPK + M6, J)+CC2*BET*KNT*BET+*2/KNTMKNT/XL < J ) )♦ (-VELF+BET* A 241 
lCC2*VELF*BPK + GC-**2 + DVX(l ) ) )/XL (J ) A 242 

K(7*J)— VELF A 243 

K(8»J)”1.0 A 244 

GO TO 18 A 245 

15 K(1#J)«-(VFLF+XL(J) ) A 247 

K ( 2* J ) *1 • 0 A 248 


VO 

ov 


K ( 3# J ) -0. 0 A 249 

K(4,J)-0.0 A 250 

K<5#J)-0,0 A 251 

K<6*J>«0.0 A 252 

K(7,J)-0.0 A 253 

K(8,J)-0.0 A 254 

GO TO 18 A 255 

C A 256 

16 BPR- (PR1*VELF + XL( J) )*XL ( J )-BPK A 257 

K<4, J) — ( VELF + 2.*Xt (J) >/(2.0*GG**?*DVX<l) ) A 25 8 

K ( 2* J ) -1 • 0 A 259 

K(l> J)-(BKN*(1«/XL( J)-1.0)+(BET*XL<J>+2.0*GG**2+DVX(1) ) ♦K f 4# J ) )/XL A 260 

1(J) A 261 

K(3> J)-<BKN/XL < J > )*K(4, J ) A 262 

K(6, J)«(GG**2+DVX(1 )-CC2*VELF*BPP)*K(4, J )/BPP A 263 

K<5> J)*-(BPK*K(6, J) + (CC2*VELF*BPK + GG**?+DVXU )+CC2*XL < J )*BFT)*K (4> A 264 

1 J ) ) / XL ( J ) A 265 

K(7>J)-(KNT-XL(J)4(VELF+XL(J))/KNT)*K(4>J) A 266 

K(8*J)*(XL(J)/KNT)*K(4jJ ) A 267 

GO TO 18 A 268 

C A 269 

17 K(1,J)-0.0 A 270 

M2#J)-0.0 A 271 

M3#J)-0.0 A 272 

K(4,J)-0.0 A 273 

K ( 5> J ) •-( PR1*VELF+XL(J ) ) A 274 

K(6,J)-1.0 A 275 

K ( 7* J ) -0 • 0 A 276 

K(8*J)-0*0 A 277 

IB CONTINUF A 278 

19 CONTINUE A 279 

C A 280 

GO TO 29 A 281 

20 CONTINUE A 282 

C************************************************** ♦*♦*** a 283 

C EIGENVECTORS FOR BFT-0 £ MOD-1 A 284 

C FOUR REPEATED EIGENVALUES A 285 

C *♦**************♦**************♦*************+♦******** A 286 



c 

WRITE ( 6# 100 ) XL1#XL2»XL3>XL4#XL5*XL6>XL7*XL8 

A 287 


DO 21 J *1 # 4 

A 286 

21 

XKJ>—KNT 

A 289 


00 22 J-5*8 

A 290 

zz 

XL ( J )«KNT 

A 291 

c 


A 292 


DO 28 J»l* 8 

A 293 


GO TO t23#24,25. 26*23*24,25*26), J 

A 294 

c 


A 295 

23 

Ml# J)«0.0 

A 29fc 


M2# JJ-0.0 

A 297 


M3# J)*1.0 

A 298 


M4#J)«0*0 

A 299 


M5# J>»0.0 

A 300 


M6# J)«0.0 

A 301 


K(7*J> — XL(J)/KNT 

A 302 


K ( 8# J)— KNT 

A 303 


GO TO 27 

A 304 

C 


A 305 

24 

M 1# J ) *1 • 0 

A 306 


K (2> J ) -XL C J ) 

A 307 


M3* JJ-1.0 

A 308 


M4# J)-0.0 

A 309 


K<5* J)«2*0 

A 310 


K ( 6# J)«2.0*XL( J) 

A 311 


M7* J)--XL< J)/KNT 

A 312 


K ( 8# J ) *— KNT 

A 313 


GO TO 27 

A 314 

C 


A 315 

25 

Ml# J)«0.0 

A 316 


K ( 2# J ) »0* 0 

A 317 


K ( 3# J ) *1 • 0 

A 318 


K(4* J)«2.0 

A 319 


M 5 # J ) *0 • 0 

A 320 


K (6# J)»0* 0 

A 321 


M7# J)«-(1.0 + XL(J))/KNT 

A 322 


M8# J)«-XM J)MXL< J) +2.0) /KNT 

A 323 


VO 

'•si 


v£> 

00 



GO TO 27 

A 

324 

c 


A 

325 

26 

K(1,J)*0*0 

A 

326 


K ( 2, J ) *0* 0 

A 

327 


K(3,J)«1.0 

A 

328 


K(4, J)«XL ( J )*( 2.0*XL ( J ) + 3.0) /KNT**2 

A 

329 


K(5, J)*4.0/GG**2 

A 

330 


K(6, J)*4.0*XL ( J )/GG**2 

A 

331 


K(7,J) — <1.0 + Xl< JM/KNT 

A 

332 


K(8, J)«-(1.0 + XL ( J) )**2/KNT 

A 

333 

27 

CONTINUE 

A 

33 A 

28 

CONTINUE 

A 

335 

29 

CONTINUE 

A 

336 


DO 30 L*1 , 8 

A 

337 


DO 30 1*1,8 

A 

338 

30 

KK(I,U*K(I,L) 

A 

339 

C 

WRITE (6* 100) < (K(I,L),L*1,8),I*1,8) 

A 

340 


CALL MATINV ( 8 , B, K , 0, C D, 0, DETR , C S , PV, IK ) 

A 

341 

C 

WRITE ( 6* 203 ) DETR 

A 

342 

C 

WRITE (6, 100) ((K(I,L),L*1,8),I*1,8) 

A 

343 


KARRAY(1)«20 

A 

344 


KARRAY(2>*8 

A 

345 


KARR AY C 3 ) *8 

A 

346 


KARR AY ( A ) *8 

A 

347 


KARR AY ( 5 ) * 8 

A 

348 


KARRAY<6)-8 

A 

349 


KARR AY (7) *8 

A 

350 


CALL MATOPS (KAPRAY,K,C,FF) 

A 

351 


CALL NATO PS ( K AP P AY, K K , K , DD ) 

A 

352 


CALL MATOPS (KARP AY, FF , KK ,CCC ) 

A 

353 

C 

WRITE (6,100) ((C(I,L),L*1,8),I*1,8) 

A 

354 

C 

WRITE (6,100) ( (DD(I,L),L«1,8),I*1,8) 

A 

355 

C 

WRITE (6,100) ( (FF(I,L),L*1,8),I*1,8) 

A 

356 

c ********* ******************************************* 

A 

357 

c 

ELEMENTS OF INVERTED MATRIX 

A 

358 

c *********** *************************************** 

A 

359 



A(1*1I-K<5,1> 

A 

360 

A{1*2)»K(5,2> 

A 

361 

A(1,3)«K( 5,3) 

A 

362 

A(l*4)-K<5,4) 

A 

363 

A(1,5)*K(5,5) 

A 

364 

A(l*6)-K(5,6> 

A 

365 

A(l,7)-K< 5,7) 

A 

366 

A<l*fi)-K<5,8) 

A 

367 

A(2*1)«K(6,1) 

A 

368 

A(2,2)-K( 6,2 ) 

A 

369 

A(2*3)-K<6,3> 

A 

370 

A(2*4)-M6,4> 

A 

371 

A(2,5)-K<6,5) 

A 

372 

A ( 2*6 ) *K ( 6, 6 ) 

A 

373 

A(2*7)-M6,7) 

A 

374 

A ( 2# 8 1 »K ( 6, 8 ) 

A 

375 

A ( 3# 1 ) -K ( 7, 1 ) 

A 

376 

A(3*2)-K(7,2) 

A 

377 

A(3*3)»K<7,3) 

A 

378 

A{3,4)-K<7,4) 

A 

379 

A(3*5)-K(7,5) 

A 

380 

A( 3*6)«K(7,6) 

A 

381 

A ( 3* 7 ) *K ( 7, 7 ) 

A 

382 

A ( 3* 8) -K ( 7, 8 ) 

A 

383 

A ( 4*1 ) -K { 8, 1 ) 

A 

304 

A(4,2)«K( 8*2) 

A 

385 

A ( 4* 3 )*K ( 8, 3 ) 

A 

386 

A ( 4* 4 ) -K ( 8, 4 ) 

A 

387 

A<4,5)«K( 8,5) 

A 

388 

A ( 4* 6 ) -K ( 0, 6 ) 

A 

389 

A ( 4* 7) «K ( 8* 7 ) 

A 

390 

A(4*8)-K(8,8) 

A 

391 

CO 31 1-1,4 

A 

392 

00 31 J-1,8 

A 

393 

B ( I * J ) ■ 0. 0 

A 

394 

IF ( INDF X • EQ • 2 ) GO TO 32 

A 

395 


vo 

vo 



C ************************** t ******** t* +t ******** t* +* A 396 

C OUTER BOUNDARY CONDITION FOR REGULAR PROBLEM A 397 

£*************************************************** A 398 

B<1,1)«1.0 A 399 

B(2,3)«1.0 A *00 

B(3,5>«1.0 A *01 

B I*, 7) »1 • 0 A *02 

GO TO 33 A *03 

32 CONTINUF A *0* 

q***************************** ********** ********** A *05 

C OUTER BOUNDARY CONDITION FOR ADJOINT PROBLEM A *06 

C***************************** ******************** A *07 

B<1,2)«1.0 A *08 

B(2,*)«1.0 A *09 

B<3,5)-1.0 A *10 

B(*, 8)*1*0 A *11 

33 CONTINUE A *12 

DO 3* I-l,* A *13 

ALPHA ( I ) ■ 0* 0 A *1* 

3* BETA(I)«0.0 A *15 

C**************************************************** A *16 

C PARAMETERS FOR JOCK A *17 

C***************************************************** A *18 

NY-IE A *19 

IFLAG-0 A *20 

AE* ACC A *21 

RE-ACC A *22 

IWORK ( 1 ) * 10 A *23 

I WORK ( 11 ) ■ 1 A *2* 

DO 35 1-1,10 A *25 

35 WORK Cl) «Y( 9*1 + 1) A *26 

XSAVF-0.0 A *27 

KL-1 A *28 

C A *29 

CALL JOCK <Z,8,8,Y,NY, A,*, AL PHA , *, 8, *, B E TA, * , 0, R E, AE , I F L AG, WORK , 10 A *30 

1000, IWORK, 150,0) A *31 
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IF ( INDEX • EQ • 2 ) GO TO 36 
MNV • IWORK ( 1 ) 

MNW»IWORK(2) 

AD«Z(2,NY) 

DA«Z(7,NY)/AD 
GO TO 37 

36 CONTINUE 
MNV* IWORK ( 1 ) 

MNW-IW0RM2) 

AO» Z f 1 » NY ) 

DA«Z(8,NY)/AD 

37 CONTINUE 

IF C ITR • NF • 1 ) GO TO 38 

CALL ITRN (GG,KNT,DA,IN,IM,ITP,IFLAG,*NV,MNW) 
GO TO 39 

38 CALL ITRN ( KNT , GG, 0 A, I N, I K f I TP , IFL AG, KNV, MNW ) 

39 CONTINUE 

IF (-IN.EQ.NIT) GO TO 60 
IF ( IN . EQ « 1 ) GO TO 3 
IF (INDEX. EC. 2) GO TO AO 
NN»2 

GO TO *1 


A 432 
A 433 
A 434 
A 435 
A 436 
A 437 
A 438 
A 439 
A 440 
A 441 

A 442 
A 443 
A 444 
A 445 
A 446 
A 447 
A 448 
A 449 
A 450 
A 451 
A 452 
A 453 


40 NN« 1 A 454 

41 CONTINUE A 455 

DO 42 J * 1 , N Y A 456 

DO 42 l«l,6 A 457 

42 Z(l,J)*Z(L,J)/Z(NN,NY) A 458 

WRITE (6,53) A 459 

DO 43 J«1,NY,5 A 460 

43 PRINT 58, Y { J ) , ( Z ( L , J ) , L «1 , 8 ) A 461 

C WRITE(6,76) A 462 

C WRITE (6, 100) { ( KK ( I ,L ),L*1,8),I*1,8) A 463 

C WRITE(6, 81 ) A 464 

C WRITE(6,100) ( (CCC (1,1 ),L*1,8),I«1, 8) A 465 

WRITE (6,59) A 466 

WRITE (6,54) (XL(I),I»1,8) A 467 

PRINT 55, IN A 468 


o 



c A 469 

NUM0RT»IW0RK(1 ) A 470 

PRINT 57, KNT, GG, BET’ A 471 

C WRITE (6,100) Z ( 1 , N Y ) A 472 

GO TO 44 A 473 

STOP A 474 

44 CONTINUE A 475 

45 CONTINUE A 476 

IF (0.2.LT.KNT.AND.0.5.GE.KNT) GO TO 46 A 477 

IF (KNT. LE. 0.2) GO TO 47 A 478 

IF (KNT. GF. 1.0) GO TO 48 A 479 

IF (0.5.LT.KNT.AND.1.0.GT.KNT) GO TO 49 A 480 

46 KNT *KNT-0 • 050 A 481 

GG-GG-0.20 A 482 

GO TO 52 A 483 

47 KNT ■ KNT— 0 .02 A 484 

GG-GG+0.1 A 485 

GO TO 52 A 486 

48 KNT-KNT-0.20 A 487 

GG-GG-0.3 A 488 

GO TO 52 A 489 

49 IF (BET. IE. 2.0) GO TO 50 A 490 

KNT •KNT-O. 05 A 491 

GO TO 51 A 492 

50 KNT-KNT-0.05 A 493 

51 GG-GG-0.3 A 494 

52 CONTINUE A 495 

IL-IL+1 A 496 

IF (IL.GE.ILM) STOP A 497 

GO TO 2 A 498 

STOP A 499 

C A 500 

53 FORMAT ( / / , 1 5HE IGEN FUNCTIONS,//) A 501 

54 FORMAT (8X,8F15.6) A 502 

55 FORMAT ( / /, 10X, 17HNQ OF I TER AT I ONS «, 1 2 , / / ) A 503 
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56 


FORMAT <8X,2HIN,10X,2HDA,12X>?HXU10X#2HX2,12X,5H$L0PE,12X*3HDEL>1 A 50* 


12X# SHI FLAG* 8X# 3HMNV> 8X>3HMNWj / ) A 505 

57 FORMAT (1 OX# AHK NT*> FI 0 . bp 5X# 3HGG», F10 . 6> 5X, AHB ET« » F10 • fc> / ) A 506 

58 FORMAT <F5.1,8E1A,5) A 507 

59 FORMAT t//»XH >8X,12HEIGEN VALUES//) A 508 

END A 509- 

SUBROUTINE I TRN (X1*X2»DA* IN> IM* ITR* IFLAG*MNV*MNW) B 1 

***************************************************************** b 2 

ITERATION BY NEWOTON RAPHSON METHOD B 3 

ACC IS THE ACCURACY 6 A 

ITR DEPENDING ON THE VALUE OF ITR ITERATION IS PERFORMED ON B 5 

KNT OR G. X 2 REMAINS CONSTANT THROUGHOUT THE ITERATION. B 6 

AL-DEL/X1 IS THE PARAMETER WHICH IS CHECKFD FOP ACCURACY B 7 

*♦ *************************************************** ************* b 8 

COMMON /E EE/ ACC>FACT B 9 

IF CIN.GE.l) GO TO 1 B 10 

DAA«DA B 11 

DEI*0. 005+X1 B 12 

X1»X1*DEL B 13 

IN-IN+l B 1* 

IM»1 B 15 

RETURN B 16 

1 SLOPE*(DA-DAA)/DEL B 17 

DEL--DA/SLOPE B 18 

AL-DEL/X1 B 19 

WRITE <6>6) IN>DA,Xl,X2iSL0PE>DEL,IFLAG*MNV,MNW B 20 

IF C ABS( AL ) .LE.ACC) GO TO A B 21 

2 IF «ABS(DEL).LE.(FACT*X1) ) GO TO 3 B 22 

DEL *OEL*C. 7 R 23 

GO TO 2 B 2 A 

3 CONTINUE B 25 

GO TO 5 B 26 

A IM*0 B 27 

RETURN B 28 

5 CONTINUE B 29 

DA A*DA B 30 


o 

to 



Xl-XUDEL B 31 

IN* IN*1 B 32 

IM-1 B 33 

RETURN B 34 

C B 35 

6 FORMAT (8X*I3>2XfE15.6f4X,F8.5#4X,F8.5>4X,E15.6»4X,F10.5>4X,I4,eX# B 36 

1 14 # 8X * 1 4 ) B 37 

END B 38- 

SUBROUTINF F MAT ( X, S# SP# IGOFX,Et EP ) C 1 

PEAL KNT»MACH#MU»MUP C 2 

DIMENSION SCI), SP ( 1 ) , E ( 1 ) » EP ( 1 ) C 3 

COMMON / AAA/ KNT>BET>VELF>GG»MACH>CC>PX>PY C 4 

COMMON /B BB / XSAVE »Kl > INDEX C 5 

C EIGHT FIRST ORDER EQUATIONS FOP REGULAR AND ADJOINT PROBLFM C 7 

C AND THEIR COEFFICIENTS C 8 

C** ***************************************************** ************ C 9 

CC1-CC+1.0 C 10 

CC 2*CC +2 • 0 C 11 

GAM* 1, 4 C 12 

EC- (GAM-1. 0)*MACH**2 C 13 

IF (XSAVE.EQ.X) GO TO 1 C 14 

IF (XSAVE.LT .X ) KL-1 C 15 

XSAVE-X C 16 

CALL PROF (X>UfUR*UPP>T>TP*TPP*PRANDL>MU*MUP*ALFA*ALFAP#DUXfDTX>DM C 17 

lUXfDALFX>V>VP>DVX>VPP*DUPX>DTPXjKl ) C 18 

PR-PRANDL C 19 

TPTX«TP*V+DTX*U C 20 

A44P*4.0*TP*TPTX/T**3-?#0*(TPP*V + DTPX*U + TP*VP + DTX*UP ) / T**2-TP* ( DUX C 21 

l+VP) / T**2M0UPX + VPP) /T + UP*BET/T-U*BET*TP/T**2 C 22 

C IN ABOVE LINE TERM -U*TP/T**2 WAS DELETED AT THE END C 23 

C C 24 

A21«KNT**2+(BET*U+DUX)/(T*MU) C 25 

A22* (V/T-MUP ) /MU C 26 

A23--UP/ ( MU*T ) C 27 

A24-0.0 C 28 


VOI 
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A25»-((U*DUX«fV*UP)/T**2 + ALFA*UPP*ALFAP*UP)/MU C 29 

A26— ALFA*UP/MU C 30 

A27-0.0 C 31 

A28-0.0 C 32 

C C 33 

A31--CBET-DTX/T) C 34 

A32«0.0 C 35 

A33-TP/T C 36 

A34-0.0 C 37 

A35«U*BET/T+<0UX+VP)/T-2.0*TPTX/T**2 C 38 

A36-V/T 39 

A37*— KNT C 40 

A38*0«0 C 41 

C C 42 

A41— DVX/T-2.0+U*GG**2/T-(DTX/T)*( V/T-CC2*MUP)+BET*(V/T-2,0*MUP)+M C 43 
1U*CC2*(TP*(DTX/T-BET)/T+PR*V*(DTX/T-FC*PX)/(MU*T)+DTPX/T-TP*DTX/T* C 44 
2*2) C 45 

A A 5* CV/T-CC2*MUP )* ( 2 . 0*T PT X / T** 2- < DUX* VP ) / T-U*BE T/ T ) +MU*CC 2* ( ( PR*V C 46 
1/(MU*T) )♦ <U*BET/T+MU*KNT**2/PR-TPTX/T**2-EC*AIFA*UP**2-TP*ALFAP/PR C 47 
2-ALFA*TPP/PR)-(TP/T)*(2.0*TPTX/T**2-(DUX+VP> /T-U*BET/T ) +A44P )+ALF A C 48 
3*<CC1*DUPX+CC2*VPP+BET*UP)MU*DVX + V*VPMU*GG)**2 )MT**2 )+DALFX*UP+ C 49 
4ALFAP*(CC*DUX+CC2*VP) C 50 

IN THE ABOVE LINE U*GG**2 WAS MODIFIED AS (U*GG)**2 C 51 

C 52 
C 53 

A48— KNT + MU C 54 

A46— V*(V/T-CC2*MUP)/T+MU*CC2*(TP*V/T**2-2.0*TPTX/T**2+(DUX+VP) /T + C 55 
1U*BET/T+PR*V*( V/T-MUP/PR-AIFA*TP/PR) / ( T*MU) +VP/T-V*TP/T**2 ) +AIF A* ( C 56 
2CC*DUX+CC2*VP> C 57 

A47»KNT*V/T-2.C*KNT*MUP-Mll*CC?*KNT*TP/T C 58 

A42-CC1*MU*9ET*0MUX+MU*CC2*( DTX / T-BE T-2 • 0*PR *V*EC*UP / T ) C 59 

A43»-U*BET/T-MU*KNT**2-VP/T-(TP/T)*( V/ T-CC2*MUP ) +MU*CC2* ( PR*V* < TP/ C 60 
1T-FC*PY)/ (MU*T)+TPP/T) C 61 

A44-0.0 C 62 

C C 63 

A61-PR*(DTX/T-EC*PX ) /MU C 64 

A62— 2,0*PR*EC*UP C 65 


o 

cn 


c 


c 

c 

1 


A63-PR*(TP/T-EC*PY) /MU 
A64-0.0 

A65-PR*(U*BET/T4-MU + KNT**2/PR-TPTy/T**2-EC*ALFA*UP**?-ALFA*TPP/PR-A 
1LFAP*TP/PR ) /MU 

A66-PRMV/T-MUP/PP-ALFA*TP/PR ) / MU 

A67-0.0 

A68-0.0 

A81*DMUX*KNT/MU*CC1*KNT*DTX/T 

A82-0.0 

A83*KNT*MUP/MU+KNT*CC 1+TP/T 
A84*-KNT/MU 

A85«ALFA*KNT*CC*(DUX + VP)/M:U“KNT*CC1*(2.0*TPTX/T**2-(DUX+VP)/T-*U*BE 
1T/T ) 

A86«KNT+CC1*V7T 
A87»U*BET/ <T*MU)+KNT+*2 
A88-V/(T*MU)-MUP/MU 

WRITE (6*100) A44R, A41,A45*A46*A61*A81*A85 
CONTINUE 

IF (INDEX. EQ. 2) GO TO 2 
SPC13-SC2) 

SP(2)«A21*S<1)4A?2*S(2)+A23*S<3 ) + A 24*S ( 4 ) ♦ A 25»$ ( 5 ) +A26*S (6 ) +A27*S C 
17 ) +A28*S ( 8 ) 

SP(3)»A31*S(1 )+A32*S(2)+A33*S(3)*A34*S(4)+A35*S(5)+A36*S(6)*A37*S( 
1 7)+A38*S ( 8 ) 

SP(4)bA41*S(1)+A42*S(2)+A43*S(3)+A44*S(4)+A45*S(5)+A46*S(6)+A47*S( 

17)*A48*S(8) 

SP ( 5 ) *S ( 6 ) 

SP(6)«A61*S(1)+A62*S(2)+A63*S(3)+A64*S(4)*A65*S(5)*A66*$<6)+A67*S( 
17) ♦A68^S ( 8 ) 

SP(7)-S(8) 

SP(8)*A81*S(1 )+A82*S ( 2)+A83*S (3 > ♦ A 84*S < 4 ) ♦ A 85* S ( 5 ) ♦ A 86* S ( 6 ) +A8 7*S < 
17)+A88*S(8) 

RETURN 


C 66 
C 67 
C 68 
C 69 
C 70 
C 71 

C 72 
C 73 
C 74 
C 75 
C 76 
C 77 
C 78 
C 79 
C 80 
C 81 
C 82 
C 83 
C 84 
C 85 
C 86 
C 87 
C 88 
C 89 
C 90 

C 91 
C 92 
C 93 
C 94 
C 95 
C 96 
C 97 
C 98 
C 99 
C 100 
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2 CONTINUE C 101 

SP(1 >— A21*$<2)-A31*S(3)-A41*S(4)-A61*S(6)-Aei*S(8) C 102 

$P(2)-— $(1)-A22*S(?)-A32*S(3)-A42*S(4)-A62*S(6)-A82*S(6) C 103 

$Pm«-A23*S(2)-A33*S(3)-A43*S< 4 )-A63*S (6 )-A83*S <8) C 104 

$P<4)— A24*$(2)-A34*S (3 ) - A 4 4*S (4)-A64*S (6)-A84*$ (8 ) C 10 5 

SP(5)--A25*S(2)-A35*S(3)-A45*S(4)-A65*S(6)-A85*S(8) C 106 

SP(6)--A26*S(2)-A36*S ( 3) -A46*S ( 4 )-A66*$ ( 6)-A86*S t 8 >-S C 5 ) C 107 

SP(7)--A27*S(2)-A37*S(3)-A47*${4)-A67*$(6)-A87*S(8) C 108 

SP(8)»-A28*S(2)-A38*S<3)-A48*S<4)-A68*S<6)-A8e*S<8)-$(7) C 109 

RETURN C 110 

C C -111 

END , C 112- 

SUBROUTINE GVEC ( X > G ) 0 1 

DIMENSION G ( 8 ) D 2 

RETURN D 3 

END D 4— 

SUBROUTINE PROF ( YARG > SU* S UP t SUP P » ST* STP* S TP P* S PR AND * SMU* SMUP* S AL F E 1 
1A*SALFP*SUX*STX*SMUX, SALFX*SV*SVP* SVX*SVPP* SUPX*STPX*KL) E 2 

REAL INTER* MU*MUP E 3 

COMMON /DDD/ Y { 101 ) » U < 101 ) * UP ( 1 01 ) * UPP ( 1 01 ) * T ( 1 0 1 ) > TP ( 1 01 ) * TPP ( 10 1 E 4 

1)*PRANDL(101)*MU(101)*MUP<101)*ALFA(101)*ALFAP<101) E 5 

COMMON /F FF / IE*MOD E 6 

COMMON /CCC/ V(101)#VP<101)*VPP(101)*DUX(101)fDUPX<101)*DVX(101)*D F 7 
1TX(101)*DTPX<101)>DMUX(101)*DALFX(101) E 8 

C***************************************************** ************ p 9 

C SUBROUTINE TO CALCULATE VALUES OF U*MU*T ETC AND THEIR E 10 

C DERIVATIVES AT A PARTICULAR Y LOCATION BY INTERPOLATION E 11 

£** ****** ****** ****** *********************** *********************** E 12 

DO 1 J-KL* IE E 13 

I-J E 14 

IF i YARG.GT.Y( J) ) GO TO 2 E 15 

IF ( YARG* EQ* Y ( J ) ) GO TO 3 E 16 

1 CONTINUE E 17 

2 MIN-I-3 P 18 

IF ( I«LE. 3 ) MIN-1 E 19 

IF ( I *GE • ( IE-2 ) ) MIN-IF-6 E 20 

SU«INTER( Y*U*YARG*6*MIN) E 21 

SUP-INTER(Y*UP*YARG*6*MIN) E 22 

SUPP-INTER (Y,UPP*YARG*6*MIN) E 23 



$T- INTER <Y,T,YARG>6>MIN) 

STP- INTER ( Y*TP»YARG>6>MIN) 
STPP-INTFR(Y,TPP,YARG>6,MIN) 

S PR AND- INTER ( Y,PRANDL, YAPG>6>MIN) 
SMU- INTER ( Y,MU*YARG#6>MIN) 

S MU P- INTER < Y,MUP*YARG>6>MIN) 
SALFA*INTER(Y>ALFA>YARG>6>MIN) 
SALFP- INTER (Y> ALFAP, YARG, 6, MIN) 

IF ( MOD • E Q • 1 ) GO TO A 
SUX- INTER (Y,DUX,YARG>6,MIN) 

SV- INTER < Y, y, YARG,6>MIN) 

SVX- INTER ( Y,DVX# YAPG*6,MIN) 

S VP -INTER ( Y, VP# YARG*6>MIN) 
SVPP-INTER (Y*VPPjYAPG*6>MIN) 
SUPX-INTER ( Y»DUPX, YARG»6»MIN) 

ST X- INTER <Y,DTX, YAPG>6>MIN) 

S MUX -INTER ( Y,DMUX, YARG#6, MIN) 
SALFX- INTER ( Y , D AL F X , Y A RG , 6 , M IN ) 

S TP X -INTER ( Y,DTPX, YAPG>6, MIN ) 

KL-I 

RETURN 

3 SU-U(I) 

SUP-UPU ) 

SUPP-UPP ( I ) 

ST-TU) 

STP-TP(I) 

STPP -TPP ( I ) 

SPR AND-PR ANOl ( I ) 

SMU-MU(I) 

S MU P -MUP ( I) 

SAIF A-ALF A ( I ) 

SALFP-ALFAP(I) 

IF (MOD.EO.l) GO TO A 
SUPX«DUPX(I) 

STX-DTX ( I ) 

STP X-DTP X ( I ) 

SMUX-DMUX ( I ) 


E 2A 
E 25 
E 26 
E 27 
E 2e 
E 29 
E 30 
E 31 
E 32 
F 33 
E 3 A 
E 35 
E 36 
F 37 
F 38 
E 39 
E AO 
F Al 
E A 2 
E A3 
E AA 
E A 5 
E A6 
E A 7 
E A 8 
E A 9 
E 50 
E 51 
E 52 
F 53 
E 5 A 
E 55 
E 56 
E 57 
E 56 
E 59 
E 60 
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SALFX-DALFXm 

SV*V<I) 

SUX-DUXU ) 

$VX"OVX(I) 

$VP»VP ( I ) 

$VPP«VPP<I ) 

KL ■ I 
RETURN 

4 CONTINUE 

SV-0.0 
SVP-0.0 
SVPP-0.0 
SVX*0.0 
SUX*0.0 
SUPX-0.0 
STX-0.0 
STPX "0 *0 
SMUX-0.0 
SALFX-0.0 
KL-I 
RETURN 
END 

PEAL FUNCTIONINTEP(X,Y,XAPG# I DEG# MIN) 
DIMENSION X ( 1 51 ) # Y { 1 5 1 ) 

1 FACTOR«I • 0 
max-min+ideg 

DO 2 J«MIN#MAX 
IF (XARG.NE.Xf J ) ) GO TO 2 
INTER-YU ) 

RETURN 

2 FACTOR -FACTOR* (XARG-X1J) ) 

YEST-0.0 

DO 4 I«MIN»MAX 
TERM»Y(I)*FACTOR/(XARG-X( I ) ) 

DO 3 J-MI N# MAX 

3 IF (I.NE.J) TEPM*TERM/(X(I)-X(J) ) 

ft YEST-TEPM+YEST 

INTER-YEST 

PETUPN 

END 


E 61 
E 62 
E 63 
E 64 
E 65 
F 66 
F 67 
E 68 
E 69 
E 70 
F 71 
E 72 
F 73 
E 74 
F 75 
E 76 
E 77 
E 78 
E 79 
E 80 
E 81 
E 82- 

G 1 
G 2 
G 3 
G 4 
G 5 
G 6 
G 7 
G 8 
G 9 
G 10 
G 11 
G 12 
G 13 
G 14 
G 15 
G 16 
G 17 
G 18- 
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